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SUMMARY 


This report discusses the work partially supported under NASA Contract NAS 
8-29662, "Segregation Effects During Solidification in Weightless Melts. " The contract 
covered the period from July 1, 1973 to June, 1975, 

During the contract period, we Unexpectedly discovered a unique segregation 
pattern and evidence of strong convection currents in GaAs crystal growth. We found, 
also imexpec^edly, some beneficial effects from vibration during GaAs growth. All 
these have clear and important implications in space processing. 

In addition, the generalized problem of determining the temperature and solute 
concentration profiles during directional solidification of binary alloys with surface 
evaporation has been mathematically formulated. Realistic initial and boundary 
conditions have been defined and a computer program has been developed and checked 
out. 

The program computes the positions of two moving (evaporation and solidification) 
boundaries, their velocities of movement, and the temperature and solute concentration 
profiles in the semi- infinite material body at selected instants of time. 

The program has the following unique features: 

• Two moving boimdaiies are involved, i. e. , the evaporative boimdary and 
freezing boundary 

• Surface evaporation, and its related effects such as material loss, evaporative 
segregation, and surface cooling due to the heat evaporation, have been con- 
sidered 

• Surface temperature is realistically determined by the combined effect of heat 
radiation, evaporative cooling, and thermal diffusion 

• Material parameters such as solid and liquid densities, specific heats, thermal 
conductivities, mass diffusivities, and latent heat of fusion or evaporation, 
can all vary with both the temperature and composition 

• Realistic phase diagrams involving curved liquidus and solidus lines are used 

Typical computed results for combined evaporation and solidification are in- 
cluded. An econotechnical model of crystal growth is presented. This model, when 
fully developed, will allow: prediction of crystal size, quality, and cost; systematic 
selecting of the best growth equipment or alloy system; optimization of growth or 
material parameters, and maximization of zero-gravity effects. 
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An important aspect of this work has been to study segregations by microprobing 
GaAs samples on transverse depth profiles from different locations on the sample 
surface. Analyses of the probe data reveal a unique segregation pattern that cannot be 
explained by existing phase diagrams. The variations of solidification velocities on 
the surface also are unexpected and could be explained by us only on the basis of strong 
convection currents in the growth melt. We therefore conclude that space-grown GaAs 
crystals by the same horizontal Bridgeman method should give more uniform and 
perfect crystals, in line with other Skylab crystal growth results. 

Another important aspect of this work has been to study vibration effects on crystal 
growth. Specifically, we have carefully and systematically documented the different 
sources of vibrations on one set of GaAs growth furnaces and, in particular, fully 
characterized the beneficial vibration effects due to an unbalanced motor. The re- 
lation of vibration to liquid mixing and space processing is discussed. The report 
also includes a section on a GaAs statistically designed experiment. This experiment 
shows that with only 14 crystals, properly planned, we can study the effects of melt 
and growth variables on solute segregation and crystal quality. We believe that such, 
or similarly designed, experiments are particularly useful in space processing where 
the sample numbers are necessarily limited and the cost of experiments is too high to 
tolerate inefficiency. 


1. INTRODUCTION 


Space processing is moving closer to reality. Bigger, better, and more uniform 
single crystals of important semiconductors and welds or brazes of improved properties 
have already been made in space, as reported in the Third Space Processing Sympo- 
sium at Marshall Space Flight Center. Although processing of structural materials may 
certainly have a profit potential in the long range, it appears that the high cost per 
pound of single-crystal electronic and optical materials makes these materials the 
most desirable contenders for immediate profitable returns from space processing. 

A selected single crystal study is, therefore, highly desirable to help us understand 
the segregation effects during solidification in weightless melts. 

Important tools for tinder standing these segregation effects are analytic solutions 
or computer programs that simulate >r predict what actually happens during space 
manufacturing. Such solutions and programs are probably necessary in space pro- 
cessing and other experiments where available time and experiment^ facilities are 
limited, the cost per sample or experiment is very high, and yet only a limited total 
number of tests or test samples can be conducted. 

Theoretical predictions often save time while computer simulation saves cost. 
Specifically, anal 5 i;ic solutions and computer programs allow us to answer many 
questions during the planning or execution of space experimeysts on material solidi- 
fication, such as learning: 

• Wliat phenomena are most important and what other phenomena are negligible 

• Which influences are favorable to our imderstanding of weightless solidifi- 
cation and which are not 

• What conditions lead to optimal combination of the favorable influences or 
elimination of the imfavorable ones 

• What sample and processing conditions should be used 

• Wliat is the best way to analyze the resultant samples for understanding a 
particular phenomenon or influence 

• How to save time and money, that is, how to maximize scientific return. 

We have developed a number of analytic solutions relating to solidification and 
evaporation (Refs. 1-3). Several important computer programs have also been de- 
veloped. Some of these solutions and programs were developed imder our Contract 
NAS 8-27891, and they are already proving useful in correlating experimental re- 
sults (Refs. 4 and 5). 

However these analytic solutions and computer programs are still in their early 
stages of development. The physical models involved are very simple and require 
considerable improvements to be used for other applications. It is, therefore, an 
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important objective of this contract to refine and improve these models and the re- 
sultant analytic solutions and computer programs. 


These refined solutions and programs are more widely useful, have greater pre- 
dictive value, and provide more accurate results. Such accuracies are absolutely 
necessary to separate the rather subtle zero-gravity effects on solidification, in the 
presence of noise due to other imavoidable or unanticipated but ever-present miscel- 
laneous effects. As a result of this continued work, more efficient space experiments 
and greater scientific returns appear possible. More meaningful solidification ex- 
periments and fuller utilization of the unique space environment may also result. 

Of course, the predicted results of our refined solutions and programs should 
first be checked with selected experiments. Another objective of this contract has 
been, therefore, to design unique experiments to correlate the numerical results to 
actual solidification processes. 


A, Review of Previous Contract ■ 

Under our NASA Contract NAS 8-27891, ’’Segregation Effects During Solidification 
in Weightless Melts” (Ref. 3), two types of melt segregation effects were studied; 
evaporative segregation, or segregation due to surface evaporation, and freezing 
segregation, or segregation due to liquid- solid phase transformation. 


These segregation effects are closely related. In fact, evaporative segregation 
always precedes freezing segregation to some degree and must often be studied prior 
to performing meaningful solidification experiments. This is particularly true since 
evaporation may cause the melt composition, at least at the critical surface regions 
or layers, to be affected manyfold, often within seconds, so that at the surface region 
or layer the melting point and other .;hermophysical properties, nucleation character- 
istics, base for undercooling, and critical velocity to avoid constitutional supercooling, 
may be completely unexpected. More comments about the importance of evaporation 
is given in Appendix A. 


To predict the segregation effects of solidification time and temperature and to 
correlate these predictions with actual experimental data, ’’normal evaporation 
equations” were developed (Refs. 1, 4-6). An evaporative congruent temperature (or 
equi-evaporative temperature) was then defined and listed for various binary or ternary 
alloys. Knowing these congruent temperatures and the solute and solvent evaporating 
rates, one can predict the type (solute depletion or enrichment) and magnitude of 
compositional or constitutional changes on the critical melt surface. One application 
of this imique temperature is to explain, predict, or plan ’’anomalous” evaporative 
or constitutional melting (on cooling) or solidification (on heating) experiments. We 
then computed for a simple model the reactive jetting forces due to surface evaporation 
and, in particular, showed that these forces can be very substantial on a differentially 
heated sample and may completely destroy the unique zero-gravity environment in 
space manufacturing (Ref. 7). In addition, these jetting forces may initiate surface 
deformation and vibration or other fluid disturbances, and may even produce some 
convection currents not normally anticipated. These studies also showed which 
sample materials are preferable, which should be avoided, and what impurities are 
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harmful in producing excessive jetting or effective as stabilizing influences. The re- 
lationship between normal evaporation and normal freezing was then considered. 

Finally, applications of evaporation to space manufacturing concerning material loss 
and dimensional control, compositional changes, evaporative purification, surface 
cooling, materials standards, and freezing data interpretation were briefly described. 

In the area of segregation due to solidification, we explained in some detail the 
normal freezing process and its successful use in the semiconductor industry. Various 
constitutional diagrams demonstrated the desirability of using nonconstant segregation 
coefficient techniques in metallurgical studies. We then stated the basic normal 
freezing differential equation, together with its solutions for cases where the liquidus 
and solidus are quadratic, cubic, high-degree polynominal, and exponential functions 
of the melt temperature. The meaning of constant segregation coefficient was dis- 
cussed, together with the associated errors due to curvatures of the liquidus and 
solidus lines and the best value of constant segregation coefficient for a given solidifi- 
cation experiment. Numerical methods for computing the normal freezing behavior 
w^re then given (see Appendix B). Finally, as an example, the steady state solidifi- 
cation of the Ni-Sn system under conditions of limited liquid diffusion and nonconstant 
segregation coefficients was described. This system was studied in the M553 ©k- 
periment on Skylab. 

B, Work Areas of Present Contract 

The present contract builds on the results of the previous contract. New results 
and improvements are presented in the following areas; 

• In the computer program; 

- Realistic surface conditions of radiation and evaporation have been in- 
corporated 

- Allow material parameters (specific heat, diffiisivities, . . . ) to vary with 
temperature and solute concentration 

- Develop a simple crystal growth model. 

• Computer results for combined evaporation and solidification are presented 

• Studies of segregation in GaAs crystals 

• A statistically designed GaAs Experiment 

• Investigation of some vibration effects on GaAs crystal growth 
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2. COMPUTER PROGRAM DEVELOPED UNDER PRESENT CONTRACT 


A. Introduction 

Under the present contract, we have extended the programs to allow for 
reasonable variation of physical properties . The approach taken is to base the values 
of physical properties upon extrapolated values of temperature and concentration, 
and then to determine the values of temperature and concentration. The process is 
then repeated by re-evaluations of the physical properties. Other modifications of 
our original program are: 1) to store the physical properties for each of the mesh 
points and to employ the appropriate quantities at each step, and 2) to recheck the 
mass and heat diffusion equations to make certain that the constancy of these proper- 
ties is not assumed. 

Another major program modification has been to include evaporation effects. 
This includes evaporation before solidification that is mathematically identical to the 
problem of simple solidification in binary alloys. After solidification starts, signifi- 
cant evaporation may still exist. We then have to deal with two moving (solid-gas 
and solid-liquid) boundaries located at y(t) for evaporation and at z(t) for solidifica- 
tion, as will be described. 

The initial and: boundary conditions have also been modified to make the 
problem more physically meaningful. One such modification is to include a surface 
heat radiative loss term involving T“^. This term affects the convergence of the 
problem and creates the need for different algorithms . As reported previously 
(Ref. 3 ) the surface cooling due to evaporation is negligible for many metallic sys- 
tems such as nickel and iron alloys, or other higher melting materials. 

To obtain solutions for realistic boimdary conditions and to include various mass 
transfer effects, numerical solutions of the partial differential equations of heat and 
mass transfer are required. We have again used the finite difference method to ob- 
tain the numerical solution. 

The boundary conditions for surface temperature include radiation cooling as 
given by the Stefan-Boltzmann equation and also include evaporative cooling for both 
components of the alloy. Raoult's law has been assumed in determining the evapora 
tion rates. At the interface it is assumed that the temperature and concentration 
relationships for each phase are given by the constitutional diagram for the alloy. 

The temperature dependence of the thermal and mass diffusion coefficients are 
allowed for each phase, 

B, Evaporative Solidification of a Binary Alloy 

Given a semi-infinite binary alloy melt, initially at concentration and 
temperature T , we consider the solidification of the alloy due to surface heat loss 
by evaporation^d radiation (Fig. 1). There are two separate regimes to be con- 
sidered. The first is concerned with temperature and concentration variations 
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Fig. 1 Evaporation-Solidification of 0. 01 Sb-Ge 
Initially at 970°C at 0.01 Second 


before solidification begins; the evaporation causes the original liquid- vapor boundary 
to change. Thus, we have a moving boundary problem. The second regime begins 
with the solidification which introduces a boundary between the freezing solid and 
remaining liquid phases whose compositions, we assume, follow the phase diagram, 
i.e. , solidus and liquidus curve relations hold. Consequently, after solidiidcation 
begins, there are two moving boundaries; one is the evaporative boundary and the 
other is the freezing or solidification boundary. 

1, Equations at the Evaporative Boundary 


We denote the evaporative boundary as x = y(t) where y(0) = 0. The evaporation 
rates in mol/m^/sec for pure solute and solvent are, respectively (Ref. 6): 


U = K 10 
e 




(M T ) 
' u s' 


A /tV 1 

V=K 10 '(M T )"2 

e 'vs' 
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where Ke = 5. 83 x 10 , Mu, My are molecular weights for solute and solvent atoms, 

Ts is the evaporating surface temperatui’e in °K, and Au, Bu, Ay, By are the 
evaporating constants for solute and solvent, respectively. If ^u and ^y are the 
solute and solvent densities, then: 


, UM C 
dy u 

dt 


VM^d - C) 


where C is the concentration at the moving boundary. 

The heat loss rate equation at the boundary due to radiation and evaporation is 
given by: 


4^ = _€ (T - UT C - VT (1 - C) 
at s u v' ' 


where € is emissivity coefficient, o' the Stefan-Boltzmann constant, and 7u and 7y 
are specific heats for solute and solvent, respectively. 




The general ambient temperature Ta may not be at the absolute zero(0 K). 

4 4 4 

The T term in the equation should, then, be replaced by(T ~ T.). 
s s « 

The equation for the rate of concentration change is: 


dc 

at 


- (U - V)C 


Since the evaporative boundary is a moving one, and since both the evaporation tem- 
perature T and solute concentration C are functions of distance x = y(t) and time t, 
i.e. , T = T(x,t) and C = C(x,t), the total derivatives obtained 

from the partials, i. e. : , 


dT 

dt 




where 


aT 

dx 


and 


ac 

dx 


are evaluated at the moving boundaries . 


Given and-^, we can Integrate for y, T, and C for the moving 


botmdary using a modified Euler method. 
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^t+At 



Ai + 1 ) 
+At 





where Is the value of the derivative at time t + A t using the value v' ^ for v. 


dX 

To determine ^ at time t and t + At requires knowledge of the distribution of 

temperatures at both times. Those at time t + A t are initially approximated by an 
extrapolation and are corrected using an approximated value of the temperature of the 
evaporating boundary with the heat diffusion difference equations . Since the change in 
temperature at the boundary is greatest due to the heat of evaporation, more itera- 
tions are applied to determine it than to the determination of temperature distribution 
by means of diffusion equations . Similar considerations hold for the determination of 


dC 

dx 


and 


dG 

dt 



I The computations of the position of the evaporation boundary [y = y(t)] , temper- 
ature (T), and solute concentration (C) at this boimdary constitute an initial value 
problem in ordinary differential equations. Thus, given yo = O, T = Tq, C = Cq* at 
time t = O, and given also the equations for velocity of movement of this boimdary 
dy/ dt, and rate of change of temperature and solute concentration, dT/dt and dC/dt, 
we can determine for selected times the values of y, T, and C. The method used is 
an iterated Euler scheme; 


'n+1 




At 


n 


(y' + y’ - ) 
m + 1 


where the initial value yn + i is taken as y^. This scheme must be connected to the 
problem of determining the temperature and solute concentration distribution within 
the semi-infinite body because the derivatives dy/dt, DT/dt, and dC/dt depend upon 
these quantities. The first step is to determine a first approximation of the temper- 
ature and solute concentration by extrapolation and then correct these values from the 
ne^ly approximated values of the boundary position and the temperature and concen- 
tration at the boundary. 

2. Start of Solidification 


To determine the time when solidification has begun, the boundary temperature 
is compared with the temperature obtained by the inverse function for the liquidus 
cuirve evaluated at the boundary concentration. K the former is greater, then solidifi- 
cation has not yet begun. If It is smaller, then solidification has begun. In order to 
avoid an exact iterative procedure to determine the Instant of solidification and to fol- 
low it up by a starting procedure for the first time interval thereafter, a simplified 
approach has been taken that introduces a small error in the evaporative boundary and 
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friiezing boundary. By allowing the temperature to be below the solidification temper- 
ature by a small amount and by assuming that the temperatures at both boundaries are 
the same, a starting value of x = z(t) of the freezing boundary is determined so that 
the loss in concentration due to solidification is compensated by the gain in concentra- 
tion at the liquidus. Given the new temperature TI2 below the temperature at which 
solidification begins, we compute CSS = FS(TI2) andCLL = FL(TI2), the correspond- 
ing solid and liquid concentrations given by the phase diagram. To determine DELZ = 
ZI2 - YI2, the distance between the evaporative boundary and solid-liquid interface, 
we assume that the solid is entirely at concentration CSS, and the liquid varies linearly 
from CLL to CC (112) , the concentration at the first mesh point x(II2) after the evapora- 
tive boundary. The total concentration is to equal the concentration in the whole re- 
gime had no solidification taken place. We ixssume it to be CL2 computed at YI? and 
to vary linearly to CC (112) at x(II2) i This jdelds the eq^uation; 

CSS * DELZ + (CLL + CC(n2)) /2 * (x(Il2) - YI2 - DELZ) 

= (CL2 + CC(n2)) /2 * (x(H2) - YI2) 

Hence: 

DELZ * [CSS - (CLL + CC(n2)/2)] = (x(II2) - YI2) * (CL2 - CLL)/2 

where, 

DELZ = (CLL - CL2)/2 * (x(H2) - YI2) ~ (CLL + CC(H2)/2 - CSS ) 

Then, 

ZI2 = YI2 + DELZ and f 

This enables us to begin the next time step with initial values for y((^), z(tg), 

% f • (^<v) = CSS. 

C^(2(tg>) “CLL. 

3. The Two-Boundary Problem- Derivative Estimation 


The equations at the freezing boundary are those given in the Grumman Final 
Report RE-458 to Contract NAS 8-27891 (Ref. 3), with the exception that the freezing 
boundary is now called x = z(t) and not x - y(t) as in Eq. 49 c-g. At every time step 
we must compute (in addition to the temperature and concentration at the evaporation 
boundary) the temperature at the freezing boundary. The concentrations are deter- 
mined by the phase diagram. The method we employ is that which determines T (the 
solidification temperature) and^ by means of Eq. (49) f,g. Having obtained^ we 
obtain z(t) by means of a modified Euler method. Since the Eq, (49) f,g required 

approximation for ^ and \^J ^ ^ we must develop techniques 


for these approximations appropriate to various situations for mesh points. In addi- 


iidx/ y,t 


and 


W/y,f 


tion, for the computation $ , and ~ at the evaporative boundary, we also need 
/dT\ / t»t et 


dt 


When there are two mesh points between y and z, then the 
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techniques alluded to above are available. This involves determining d^T at both y 
! 2 

apd z and the same for 9 C . When there is only one mesh point between y and z, then 
2 

8 T at both points are the same. When there are no mesh points between y and z, then 

_ T(z) - T(y) 


dx 2 

we can assume either that 9 T is zero and hence 

dx 


(» 1 ) 

Vdx/y \dx /i 




1 .+ 


Z - y 
z - y 


or that 


dX rt n 

A^'v d"T 2 dT 

choice of k must be small so that -X— k = k - — is negligible. Thus/ since 

a/2 ax 

9^T 

z-y is also very small this option is indistinguishable from -—r = 0. We have three 


2 


:=o( 


£T\ 

dx/y 


9 ^ d\ 

qases: 1) no mesh points between two boundaries and we assume 

o dx^ 

( dT \ d^T 

-—I ; 2) one mesh point between y and z when — ^ is obtained from the three 
9z /z _ o dx ^ 


points and 


■ (m- 


T(z) - T(y) (z ^ y) d^T 


z-y 


dx 2 


and 


/dT\ 

\d X /z 


T (z) - T(y> 
z-y 


(z - v) ^ r, and 3) when two or more mesh points, say x. and x. are between 
•V and z so that we can computef —V and( ^ T | separately and distinct. Then: 


/_ar\ 
V dx/j 


/ d^T\ 

( “■'■■o ) . In general, it is necessary to compute- - ^ and in three ways; two 
\ dx‘“/z dx^ 

ways indicated above for the solid regime and a third for the liquid side of the freez- 
ing boundary. It is similarly necessary to compute_|!Z.and.^.C in three ways. 

4. Boundary and Mesh Points 


_ T(Xi) - T(y) _ 

Xj - y 


(x,--y) / d^TV , T(z)-T^\+ l 

2 \9x^ ^-^i + 1 


“ ^i + l) 




d^C . 


When boundary points come close to mesh points, the computation of derivatives 
may be vitiated by closeness to mesh point. Therefore, tests are made to determine 
when such closeness occurs as usually expressed in terms of a decimal fraction of 
the interval. In that case; the reference point is moved to the next mesh point and the 
values of T and C at the skipped mesh point are obtained by linear interpolation. This 
interpolation depends on which side of the solid-liquid Interface the mesh jx)int lies . 
Por the evaporative boundary similar considerations hold. 


5. Solution for Remaining Points 

The solution for the remaining ix)ints is obtained as in the Final Report pre- 
viously mentioned, pages 3-14 and 3-15 (Ref. 3). One change is, however, necessary 
because the first mesh point (or more) are no longer under consideration if the evap- 
orative boundary has passed them. The subroutine TRKT is used to solve for the 
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remaining points . In this subroutine we compute the values of temperature and con^ 
centration at intermediate mesh points when given the values at the two extreme mesh 
points. We replace the values at the mesh point to the left of the evaporative boundary 
by those at the evaporative boundary point, before solving for the intermediate points. 
This can be done without destroying any useful information since that mesh point is no 
longer used in the computations. The subroutine TRET does not depend upon equal 
spacing or any regular spacing and therefore can accommodate this usage. 

6. Convergence 

The convergence problem is the crux of the program. Oscillation tends to 
cause the needed quantities to overflow. Thus, tests must be made on all the quan- 
tities to contain them within reasonable boimds. The subroutine MOTON is used to 
check the monotonicity of these consecutive points. In addition, the temperature at 
the evaporative boundary is necessarily less than the temperature at the freezing 
boundary. This condition is always imposed in the program. 

In addition, the solution for the solidification temperature and freezing boundary 
derivative (especially the latter) involves very rapidly changing quantities . More 
iterations should, therefore, be expended in this part of the program. Fewer itera- 
tions are needed for determining the evaporative boundary, and the temperature and 
concentration at that boundary. The program allows five Iterations in the former for 
each of the latter. The number of iterations of the latter is used in a manner anal- 
ogous to that described in Ref. 3. 

An input quantity NIT (usually a multiple of 4) gives the maximum number of 
iterations. When NIT/2 iterations occur and convergence is not reached, the time 
step size is halved. This process is continued until either convergence is attained 
or the minimum step allowed by the program has been iterated NIT + 1 times. In this 
case the program may stop or continue on using the nonconverged quantities. Very 
often these quantities are sufficiently smooth so that convergence will occur on the 
next interval and the program gives satisfactory results . 

However, if the program proceeds with the minimum step and the maxtinuTn 
number of iterations, the results may be spurious. In case of overflow, there is no 
doubt of it. Otherwise the user must look at results to decide whether he finds them 
reasonable. 

7. Eutectic Temperature 

When the Interface temperature T| reaches the eutectic temperature T^, the 
program is changed so that interface temperature from then on remains at T^. The 
rate of change of the boundary position is etill determined from Eq. 49g on page 16 
of Ref. 3. This leads to a rapid growth in the interface boundary. 

C. Variable Parameter Program Modification 

The original evaporative solidification program as well as its predecessors 
assumed that all the physical parameters were constant. Recognizing that this assump- 
tion is often a poor model of the real situation, we introduced variable parameters as 


functions of melt temperature, or solute concentration, or both. As an Initial effort 
in this direction, we limited the arbitrary functions to linear ones in the absence of 
any better known functional forms , The basic program logic for the linear approach 
would not change if non-linear terms are Introduced in the future. 


The specific heat of the solute or of the solvent Cy are assumed to be linear 
functions of temperatures T only (Ref. 8): 


C ~ C + C *T and C = C + C *T 
u uo uT V VO vT 


The density of the solute, p and of the solvent, are also assumed to be linear 


functions of temperature above; 


V 


p = p +P ^ *T and P = p +p ^»T 
u uo uT V VO vT 

The mass diffusion coefficients for the solid D is assumed to be a linear func- 
tion of both concentration and temperature: 

D =D +D *C +Dm *T» 
s so sc s sT 


The mass diffusion coefficient for the liquid is similarly defined: 

The thermal conductivity of the solid ksu» kgy and of the liquid k/ ^ k^v for solute and 
solvent, respectively, are assumed to be linear functions of temperature: 

k =k , +k t,«T, k =k +k rr*T, kp =k/ +kp „.T, k/ =kp +kp r^rT. 

sU suo suT ’ sv svo svT * Ju -*uo >«uT ’ xyo -^uT 


The thermal diffusion coefficients a ^ and Ug are determined from the equation , 
a^^ = ky /p c and = ks/p c, where k^ = c^ k^u + (l-c^)'k^v and kg = Cgkgu + 

(1-Cg) k^^, and Pand c are obtained from p p y and c^ and Cy respectively. Thus 
P = Pu c + p y (1-c) using the value of concentrati'^n for liquid, or Cg for solid. 
Similarly c = c,-,^ C + c (1-C) (Ref. 9) using the value of concentration Cy for liquid or 
Cg for solid. 

With these changes in physical parameters certain changes in program follow. 

In the original programs a subroutine TRIST was used to determine the concentration 
or temperature distribution at the new time t + At, given the values at time t and the 
values for the extreme points at t + A t. This assumed a constant diffusion coefficient 
on either side of solid-liquid interface. With a variable diffusion coefficient the pro- 
gramming gains by the use of a subroutine are lost and the four uses were converted 
to four separate sets of instructions within the main program. 

The effect of including functions of concentration and temperature Instead of 
constants make the problem of convergence more difficult. Thus more iterations or 
smaller time intervals are required to bring out convergent results. 
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D. Single Crystal Growth Prc^ram 


The computer program for crystal growth is a variant of the non-evaporative 
solidification program with a number of basic modifications j histead of the three 
options for surface temperature, only the linear surface temperature, is provided. 

The conditions at the solid-liquid interface are changed by assuming only the equation 
for the redistribution of temperature at the moving boundary and replacing the one on 
concentration by assuming a constant derivative for the moving boundary. This sim- 
plifies the interface problem because it leaves only one imknown, the interface tem- 
perature. The two conditions of linear surface temperature and linear varying inter- 
face boundary simulate the conditions of crystal growth, the constant temperature 
gradient about the linear varying interface boundary. 

The single crystal growth model is shown in Fig, 2. The computer program 
for this model assumes a constant temperature gradient and a constant (solid-liquid) 
interface velocity. This program is a variant of the linear surface temperature pro- 
gram developed in 1972 imder our contract NAS 8-27891 (Ref. 3). In this program, 
the constant (solid-liquid) interface boundary velocity replaces one of the equations 
involving the interface boundary derivatives. This simplifies the basic iteration 
problem because now only the interface temperature must be obtained in a manner 
consistent with heat and mass diffusion equations. The linear surface temperature 
propagates a nearly linear temperature profile at the interface which, combined with 
the constant growth velocity, gives a nearly constant temperature gradient about the 
solid-liquid interface boundaiy. 

Tables 1 and 2 indicate the insensitivity of the evaporation boundary and its 
temperature to grid spacing, provided that the spacing is always larger than the 
evaporation boundary point. Tables 3 and 4 Involve solidification boundary and show 
that the temperature is insensitive to DELX but that the solidification boundary is 
quite sensitive to the choice of DELX. Thus, it is important to use DELX sufficiently 
small so that the solidification boundaiy movement is fully exhibited and not stunted 
by a large grid spacing relative to which the boundary size is small. The spacing 
affects the evaluation of the first and second temperature partial derivatives with 
time which are larger in absolute values for sm^ler spacings , due to more rapid 
temperature changes near the boundaries . 

The figures (Figs. 3-5), prepared from the computed results in Tables 1-4, 
indicate the superiority of unequal over equal grid spacing. For DELX = 0.01 cm, 
where the spacing is coarse, little difference is found in the temperature distribution. 
For DELX = 0.001 cm, there is greater difference between the two because the equal 
spacing has limited the semi-infinite body to 0.027 cm and fixes the temperature 
at the end point to 970°C, thus not aEowing the temperature to decline as rapidly as 
it should. For DELX = 0.0001 cm, the equal spacing method couldnot work at all 
because 0. 0027 cm is too small a range to define a semi-infinite body even for 
the small time constants under consideration. 

Fig. 4 (grid spacing DELX = 0.001), shows wide disparity between equal and 
unequal spacing whereas Fig. 3 (DELX = 0. 01) shows good agreement. The 
smaller DELX needs more points to simulate the semi-infinite, one dimensional 
case and when restricted to IX = 28, fails to allow temperature away from the 
evaporating surface to decline rapidly because it is artificially pegged at x = 0.028 
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1. Initial conditions; 

liquid cone = c^ 
temperature = 

2. Growth conditions; 

Constant growth velocity = YD 
Constant temperature gradient = G 
Constant hot side temperature = T_ 

Constant cold side temperature = 

Constant length of intermediate temperature of gradient = ZL 
Curved liquidus and solidus 

3. Wanted; 

Concentration profiles, Ge or Si crystals. 



Fig. 2 Single Crystal Growth Model 



Table 1. Variation of Temperature ( C) at Evaporative Boundary 
for Various Computation Schemes 


Scheme 

I 

n 

] 

n 

rv 

V 

Grid 

0.01 

0.01 

0. 

001 

0.001 

0. 0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 

i 

i "■ " 





0.2 

i 

966.5 

966.5 

966.5 

966.5 

966.5 

0.6 

959.4 

|959.4 

959.4 

959.4 

959.5 

1.4 

945.7 

‘ V' 

945.7 

945.7 

945.7 

945.8 

1.8 

938.9 

938.9 

938.9 

938.9 

939.0 

2.0 

935.5 

i 935.5 

935.5 

935.5 

935.6 

2.05 

934.7 

934. 7 

934.7 

934. 7 

N.C. 

2.075 

934.2 

934.2 

934.2 

934.2 

N.C. 

2.0875 

934.0 

934.0 

934.0 

934.0 

N.C. 

2.09375 

933.9 

933.9 

933.9 

933.9 

N.C. 

2.1 

933.8 

933.8 

933.8 

933.8 

933.9 

2.1125 

933.6 

933.6 

933.6 

933.6 

933.7* 

2.1375 

933.2 

933.2 

933.2 

933.2 

933.3* 

2.1875 

932.4 

932.4 

932.4 

932.4* ; 

932.5* 

2.2875 

930.8 

930.8 

930.8 

930.8* 

930.9* 

2.4875 

927.7 

927.7 

927.2 

927.7* 

927. 7* 

2.8875 

921.5 

921.5 

921.5 

921.5* 

921.5* 

3.6875 

909.3 

909.3 

909.3 

909.3* • 

909.3* 

5.2875 

885.8 

885.8 

885.8 




’*‘Hahd interpolations 
N. C / not computed 


















Table 2. Variation of Position (mhi) of Evaporative Boundary 
With Time for Various Computation Schemes 


i 


Scheme 

I 

n 

m 

IV 

V 

Grid 

Spacing 

O.Ol 

equal 

0.01 

imequal 

0. 001 
equal 

0. 001 
unequal 

0.0001 cm 
miequal 

time, ms 



: 

i 


0p2 

0.122 

0.122 

0.122 

0.122 

0.122 

0.6 

0.351 

0.351 

0. 351 

0.351 

0.351 

1.4 

0. 752 

0.752 

0.752 

0. 752 

0.752 

1.8 

0. 927 

0.927 

0.927 

0. 927 

0.928 

2.0 

1.009 

1.009 

1.009 

1.009 

1.010 

2.05 

1.029 

1.029 

1. 029 

1.029 

N. C. 

2.075 

1.039 

1.039 

1.039 

1.039 

N.C.i 

2. 0875 

1.044 

1.044 

1.044 

1.044 

N.C, 

2.09375 

1.046 

1.046 

1.047 

1.047 

N.C. 

2.1 

1.047 

1.047 

1.047 

1.047 

1.050 

2. 1875 

1.048 

1.048 

1.048 

1.048* 

1.052 

2. 2875 

1.049 

1.049 

1.049 

1.048* 

1.053 

2.4875 

1.051 

1.051 

1.051 

1.050* 

1. 055 

2. 8875 

1.055 

1.055 

1.055 

1.054* 

1.059 

3. 6875 

1.061 

1.061 

1.062 

1.061 

1.066 

5.2875 

1.072 

1.072 

1.072 




Hand interpolations 
N. C. not computed 
















Table 3. Variation of Position (Mm) of Solid-Liquid Interface 


Scheme 

I 

n 

m 

IV 

' ■ V 

Grid 

Spacing 

0.01 

equal 

0.01 

unequal 

0.001 

equal 

0.001 

unequal 

0.0001 cm 
unequal 

time, ms 






0.21 

0.205 

0.204 

0.109 

0.109 

0.106 

0.24875 

0.211 

0.211 

0.179 

0.180 

0.408 

0.28875 

0. 229 

0. 229 

0.350 

0.350* 

1.03 

0.36875 

0.283 

0.283 

0.866 

0. 860* 

2. 76 

0.52875 

0.429 

0.429 

2.291 

N.C. 

N.C. 


Table 4. Variation of Temperature (°C) at Solid-Liquid Interface 


Scheme 

I 

n 

m 

IV 

V 

Grid 

Spacing 

0.01 

equal 

0.01 

unequal 

0.001 

equal 

0. 001 
unequal 

0. 0001 cm 
unequal 

time, ms 






0.21 

933.8 

933.8 

933.8 

933.8 

t ' 

933.9 

0.24875 

927. 7 

^ 927. 7 

927.7 

927.7* 

927.9* 

0.28875 

924.5 

' 921.4 

921.5 

921.4* 

921.6* 

0.36875 

909.3 

909.3 

909.3 

909.3* 

909.4* 

0.52875 

885.8 

885.8 

885.8 



* Hand interpolations 
N.C. not computed 
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Melt Temp 




Position from Original Melt Surface, cm 


Fig. 3 Computed Temperature Distribution in Semi-mflnite Body at 
Different Times, t (initial grid spacing = 0.01 cm) 



Position from Original Melt Surface, cm 

Fig. 4 Computed Temperature Distribution in the Seml-Ihfinlte Body at Different 
Times, t (initial grid spacing = 0. 001 cm) 
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960 • 
950- 
940 - 
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I* 
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880- 
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0 


Initial Melt Temp. = 970°C 
Time t = 0.002 sec. 


t = 0. 00215 sec. 


t == 0. 0040 sec. 


» I 1 1 1 I I 

0.0002 0.0004 O.OOOo 0.0008 0.0010 0.0012 0.0014 0.0016 


Position from Original Melt Surface, cm 

F^. 5 Computed T^perature Distribution In Semi-Infinite Body at 
Different Times, t (initial grid spacing == 0. 0001 cm) 


Cm to 970°C, The unequal spacing needs but six points to give equivalent exten- 
sion and, when given 10 or 11 points, can adequately span a sufficient distance to 
sinaulate semi-infinity. At smaller DELX (0.0001) one cannot even attempt to use 
equal spacing without modif 3 ring the behavior at the last mesh point. For unequal 
spacing^ 16 ]:K>ints will adequately represent the semi-infinite body for the times 
under consideration. 

E« Programming Details 

Some details of the programming techniques and a developed program are given 
in Appendices C and D. i 
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3. COMBINED EVAPORATION AND SOLIDIFICATION 
TYPICAL COMPUTED RESULTS 


We previously pointed out (Ref. 3) that evaporation always precedes freezing 
segregation and must be studied before meaningful performance and analysis of solid- 
ification experiments. This is particularly true since evaporation in the critical sur- 
face regions may cause the melt composition to be drastically altered within times as 
short as milliseconds. Thus the melting point and other thermophysical properties, 
nucleation characteristics, and critical velocity to avoid constitutional supercooling 
of the surface region may be completely unexpected (Ref. 3). 

Evaporation, for example, may even introduce "anomalous" evaporative or 
constitutional melting (on cooling) or solidification (on heating) in simple binary 
melts. The final report to the previous contract completely discussed suah newly 
predicted phenomena. The same report also summarized and presented some simple 
analytic equations for predicting the kinetics of segregation due to surface evapora- 
tion. 

However, the study of combined evaporation and solidification was not under- 
taken under the previous contract. Nor, to our knowledge, has it been done by any 
other worker in the field in a systematic manner. Therefore, it was an important 
objective of this contract to study such a combined phenomenon. 

Specifically, we have set up the necessary differential equations to include the 
phenomenon of evaporation in our previously developed, unidirectional solidification 
model. New computer logics have been devised and computer programs prepared 
and debugged. These work areas have been described in a previous section of this 
report. In this section, we present only some typical computed results on combined 
evaporation and solidification. 

A. A Simple Example 

For our first example we computed results for the semi-infinite Ge-Sb system 
having an eutectic temperature of 590°. These results, given in Table 5, show that 
for a 0. 10% alloy of Sb in Ge initially at 970°C, directional solidification, from sur- 
face cooling through radiation and evaporation, begins shortly before 2.1 mS, with 
an interfacial temperature of 934°C and with the evaporation boundary also at 934°C, 
but located at only 0. 105 /<m from the initial melt surface. At first, the evaporating 
surface is losing material rapidly at the rate of over 50 A»m/sec, due mostly to the 
loss of Sb. After a few milliseconds, the solid surface is practically all Ge, with a 
surface receding rate of less than 1 Miu/sec, After some 30 ms, the evaporating 
boundary is practically stationary. These results agree with rough computations 
based on the evaporating rates of Ge and Sb at 700-940°C. 

The rate of solidification increases with time. Initially, the rate is very low, 
less than 0. 054 Atm/ sec on the average in the first 2. 1 ms. This growth rate in- 
creases from 2 to 12 ms due to increasing temperature gradient in the solid. A 
temporary maximum growth rate of over 0. 5 A»ni/sec is reached at 12 ms, coinciding 
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Table 5 , Evi^ration and Solidification of 0. 1% Sb in Ge Alloy Initially at 970°C 


Time, t, 
ms 

Position of 
Evaporation 
Boundary 

Temp at 
Evaporation 
Boundary 
Ty,°C 

Position of 

Solid 

Boundary 
z. A* m 

Aver 

Growth 

Rate 

Mm/sec 

Temp at 
Solid 
Boundary 
Ti, OC 

Note 

2,1 

0. 1050 

933.9 

0. 109 

54 

933. 9 

Beginning of 
Solidification 

4,0 

0.1067 

904.7 

1.96 

490 

904. 8 


8.0 

0. 1087 

847.2 

19. 19 

2,400 

849.9 


12 

0.1096 

794.3 

39. 79 

1,990 

806.6 


16 

0.1099 

748.1 

48.53 

3, 030 

750.4 


20 

0.1101 

698.9 

49.92 

2,500 

704.9 


24 

0.1101 

655.2 

51.35 

2,140 

661.0 


'28 

0.1102 

613. 5 

52.83 

1,890 

619.5 


32 

0.1102 

573.5 

65. 24 

2,040 

590. 0 

Note Large Change 
in Solid Boundary 

36 

0.1102 

534.9 

159.5 

4,430 

590.0 


40 

0. 1102 

497.5 

272.2 

6,810 

590.0 





a maximum solid temperature gradient. A temporary minimimi gro\vth rate is 
reached shortly before reaching eutectic temperature of 590°C when the solid-liquid 
interface is sufficiently inside the origin melt surface, and, hence, the effect of heat 
loss througli surface radiation is relatively small. Right after reaching the eutectic 
solidification stage, the growth rate increases rapidly because there is then no longer 
the mechanism of heat loss through sensible heat changes (since the temperature of 
the eutectic solidification is fixed at 590 C). For still longer solidification periods, 
we expect the growth rate to decrease again, because of the leveling of the solid tem- 
perature gradient. (Part of the results are shown in Table 5). 

B, Effect of Melt Composition 

This example show^s clearly the complicated process of solidification in binary 
alloys, and also the intricate interactions between surface evaporation, internal • 
solidification and heat or mass transfer. Simple solidification theories and approxi- 
mate solidification equations may, therefore, not be accurate in evaluating the cri- 
tical initial phase of nucleation and freezing. Specifically, predicting whether con- 
stitutional supercooling occurs or whether single crystals will grow may require such 
numerical computer programs as we have developed under this contract. 

To study the effect of melt composition on combined evaporation-solidification, 
we made comparative runs for the same Ge-Sb system at the same initial temperature 
of 970®C, but with the Sb concentrations originally at 0.01, 0.001, and 0. OOOl^t* At 
short times, t ~2, 1 ms, the temperature at the evaporating boundary for the four 
initial Sb concentrations is fairly constant at 934. 0 ± 1. 3°C. The position of the 
evaporation boundary, y, however, decreases proportionately with decreasing Sb con- 
centration except at the highest Sb concentration. For example, decreasing Sb from 

0. 1 to 0. 01% reduces y by a factor of 20 (instead of 10) at 2. 1 ms. After about 26 ms, 
the evaporation boundary is stabilized in position. The percentage difference in y 
between 2. 1 ms and 26 ms decreases logaritlmilcally with Sb concentration, from 5% 
at 0. 1% Sb to 18% at 0. 0001% Sb. These results again confirm the fact that evapora- 
tion of Ge-Sb alloys is mainly due to the loss of Sb. 

All these Ge-Sb alloys, according to the computed results, begin to solidify 
eutectically mthin 32 ms. Except for the 0. 1% alloy, the interfacial temperature at 
any given time is independent of the initial Sb concentration. The 0. 1% alloy has con- 
sistently higher interfacial temperatures than the more dilute alloys, by as much as 
17°C at intermediate times. The position of the solidification boundary at a given 
time is consistently smaller for the 0.1% alloy than for the more dilute alloys. The 
difference in this position increases, as expected, Avith the time. 

These and other results on combined evaporation-solidification have given us 
some unexpected data that can be summarized as follows: 

1. Surface evaporation, internal solidification, and mass or heat transfer 
interact in complicated manners. Individual treatment of only one aspect 
may not be sufficient. 

2. The rate of evaporation is very appreciable at first, but rapidly (in milli- 
seconds) becomes negligible. 
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3, Evaporation leaves permanent and extremely large segregation effects that 
can radically affect melt undercooling and nucleation. The distance of the 
evaporation controlled layer is, however, too small (at least for Ge-Sb 
alloys) to be detected by the electron probe. 

4. The rate of solidification does not monotonically increase or decrease, as 
predicted by simple, existing solidification theories; it fluctuates up and 
down. 

C. Combined Effects of Liquid Mass Transfer and Melt Composition 

The segregation effects during solidification in weightless melts are at least 
partly due to changing effective liquid mass transfer coefficient (i. e. , solute trans- 
port rate across unit concentration gradient in the liquid), I^. Various Skylab ex- 
periments have shown to decrease by up to 50 times over ground control experi- 
ments, (Ref. 2), depending on the sample shape, size, temperature, composition, 
type of experiment, and the like. Our computer program allows us to readily sim- 
ulate the evaporation- solidification behavior during the single crystal growth under 
varying D . The computed results may, therefore, even give us some insight into 
the effects of sample shape, size, temperature, composition, and type of experiments 
on weightless solidification. 

Hence, we have made nine systematic computer runs by simultaneously varying 
the initial antimony concentration c^ in germanium single crystal over three decades 
from 0.01 to 0.0001 and also owr three decades from 10"'^ to 10”® cmvsec). 

The computed results are given in Tables 6 through 9. 

As shown in Table 6, the solidification velocily z is predominantly determined 
by . Decreasing by one or two orders of magnitude reduces z by the same 
amount. In fact, the solidification velocity z can be approximated by z «» 2000 
during the time interval 0. 002 to 0. 01 second. Hence, tne solidification velocity in 
Skylab samples probably was much smaller, by up to 50 times also, during some of 
the solidification processes. Table 6 further shows that during the initial and final 
stages of solidification, the solidification velocities are further reduced by mechan- 
isms other than those related to . Specifically, in the initial stages of solidifica- 
tion, (for time t ^ 0. 001 second), the velocity is reduced about 40 times at initial 
concentration Cq = 0. 01 for all values of D/ and at higher concentrations by about 35 
times at = 10“® cm^/sec. During the later stages of solidification (t > 0,02 sec- 
ond), the solidification velocity is reduced by about twice at D^ = 10“'^ cmVsec. V 

Table 6 shows that the receding rate of the evaporating (gas-solid) boundary y 
decreases constantly with the evaporation-solidification time t, more so at high initial 
antimony concentrations c© than at low concentrations. The decreases in y, for ex- 
ample, between 0. 001 and 0.01 second are, respectively, 160, 6, and 6 times at 
Co = 0. 01, 0, 001, and 0, 0001. A decrease iu, D^ also reduces y, by about one order 
of magnitude when Cq is reduced from 0. 01 to 0, 001 during the early stages of the 
evaporation- solidification process (t L 0.002 secohd), but with no decrease when Cq 
is further reduced from 0. 001 to 0, 0001, or when t exceeds 0.002 second. Thus, the 
ev^)oration behaviors of Skylab samples were also probably maT;^edly different from 
those of ground samples. '.fi- 

A summary of some findings on the position of the solidification boundary* z, 
(Table 6) is given below. 
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c 

o 


0.01 

0. 01 

m 

2, 

cm /sec 

lO""^ 

10-5 


t, sec 



z 

0.001 

5. 0 X lO'^ 

5. 0 X 10“^ 

cm/sec 

0.002 

1.5 X lO"^ 

1.6 X 10“^ 


0.005 

2.1X10"^ 

1.7x10“^ 


0.01 

1.8 X 10“^ 

1. 8 X lO"^ 


0.02 

7.8 X 10“^ 

1.8 X 10"^ 


y 

cm/sec 


cm / sec 


t, sec 
0.001 
0.002 
0.005 
0.01 
0.02 


5. 0 X 10 


1.7 X 10 


9. 2 X 10 


3.2 X 10 


3,1 X 10 


5. 0 X 10 


1, 8 X 10 


9. 6 X 10 


3.2x 10 


Table 6. Variation of z With Dj , c t and 
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0.01 

0.001 

0. 001 

10"® 

rt* 

1 

O 

lo’® 

5.|0 X lO"® 

l.TxlO"^ 

1.8x10"® 

2. 4 X 10"® 

1.8 xio"^ 

1.8x10’® 

2. 0 X 10“® 

2.7 X 10’^ 

1.8x10’® 

1.9x10’® 

1.9xl0~^ 

1.8x10’® 

1. 9 X 10"® 

7.8 xlO’^ 

1. 8 X lO’® 
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5. 0 X 10 
1.8 X 10 

1.0 x io" 

3. 6 X 10‘ 
3. 7 X 10‘ 


8. 5 X 10 


6.9x 10 


3.6x 10 


1*2x10 

1,2x10" 


0.001 


lo"^ 


8. 5 X 10‘ 
6. 7 X 10“ 
3. 6 X lO" 
1*2X10“ 
1.2x10“ 
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01 

0.001 

0.001 

0.001 

0.0001 

o.opoi 

oi.oooi 

CO 

1 _ 

10"“* 

10 -® 

10-« 

H* 

o 

1 

4^ 

10-5 

10 “® 

X 10“^ 

1. 7 X 10"^ 

1.8 X lO"^ 

5.7x10"^ 

' -1 
1.8x10 

1.9 xio”^: 

/ 

/> 

Is- 

6.2x10? 

X 10“^ 

1.8 X lO"^ 

1.8,xl0"^ 

2.5 X ip"^ 

l.Sxlo’^ 

1.8 x io"^ 

2J5 X 10”^ 

X 10"^ 

2.7 X 10"^ 

1.8x10"^ 

2.1 X 10"^ 

2.8x10"^ 

1.8x10’^ 

2.1x10“^ 

X 10“^ 

1.9xl0“^ 

1.8 X 10"^ 

2. 0 X 10’^ 

2. 1 X 10"^ 

. 

1.8x10"^ 

2. 0 X 10"^ 

X lO"^ 

7. 8 X 10“^ 

1. 8 X lO’^ 

.. — „ 

1. 9 X 10"^ 

7. 9 X 10"^ 

1.8x10“^ 

1. 9 X 10“^ 
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mM 

0. 001 

0.001 

0.001 

0.0001 

0. 0001 

0.1)001 

Si 

-4 

10 ^ 

10"*^ 

10-5 

10"^ 

10-5 

10-« 

> 

> -4 

X 10 ^ 

8.5xi0”P ; 

^ 8. 5 X 10"® 

8.5x10"® 

8.5xl0"® 

8.5xl0"® 

8. 5 X 10"® 

;x 10"® 

6.9x10"® ; 

6.7x10"® 

7.1x10"® 

6.9x10"® 

6.7xl0"® 

7.1x^0"® 

X 10 ^ 

3.6x10"® 

3.6 X 10"® i 

3.9 xio"® 

3.6x10"® 

3. 6 X 10"® 

3.9xl0"® 

;x io“® 

1. 2 X 10"® 

1. 2 X 10"® 

l.]4xl0"® 

1.2x10"® 

1.2x10"® 

1.4 X 10"® 

ix lo”'^ 

1, 2 X 10"“^ 

1.2x lO""^ 

1 -7 

1.4 X 10 

1.2x10"'^ 

1. 2 X 10"*^ 

1.4X10"’^ 
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Table 8. Variation of z (cm) With , t, and c 


c 

o 

0.01 

0. 01 

0.01 


0.001 

0.001 


1 

o 

10-® 

10-« 


10-^ 

10 -« 

t, sec 







0.001 

^6 

1.1x10 ° 

1. 1 X lo"^ 

1.1 xlo“® 

l.lxio"^ 

1.3x10"® 

4. 8 X lO"*^ 

0.002 

1.4 X 10"^ 

1.7xl0"® 

3. 7 X lO"® 

1. 8 X 10“'^ 

2.1 X 10"® 

3.7x10"® 

0.005 

6. 8 X 10”^ 

6.6 X 10"^ 

1.0 XTO’^ 

7. 9 X 10"^ 

7.3x10”® 

1.0x10"® 

0.01 

1.6xl0"^ 

1. 5 c lO"^ 

2.0x10"^ 

1.6 xlO"^ 

1.6x10"^ 

2.0x10"® 

0. 02 

2. 6 X 10"^ 

3. 3 X lO"^ 

3. 9 x lO"'^ 

2. 6 X 10"^ 

3.4 x lO"'^ 

4.0 X 10"® 


} 

I 

I 

Table 9, Effect of c^, D^, and t on the Position of Evaporation l! 
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i 


c 

o 

0.01 

0.01 

0.01 

r 

0. 001 

0.001 

0.001 

1 

Dj 

10"^ 

10-5 

10-« 

10’^ 

10-5 

W-® 

1 

t, sec 



■ 




■3 

j 

- j 

0. 001 

5.6 X lO"*^ 

5.6 X 10 ' 

5.6x10”'^ 

5.1x10"® 

5. 1 X lO"® 

5. 1 X lO”® 

J 

1 

0.002 

5.8x10"'^ I 

5. 8 X lo"'^ 

5.8x10^'^ 

5.8x10"® 

5. 9 X lO"® 

5.9xl0"® 

i 

0.005 

' -7 ' i 

6.2x10 

6. 2 X lO"'^ 

6.2 x; 10^"^ 

7.4 x lo"® 

7.4x10"® 

7.5x10 ° 1 

,| 

0.01 

6. 5 X lO”*^ 

6.5 X lo”"^ 

6.5x10”'^ 

8.4x10"® 

8. 5 X 10“® 

.! 

8. 7 X lO”® ! 

■■ \ 

j 

0.02 

6.6 x lO"^ 

6.6 X lO"^ 

6.7x10"'^ 

8.9 X 10"® 

9.0x10"® 

9.2x10"® 

j 




































Table 8. Variation of z (cm3 With , t, and 


i 

0. 001 

n 

0. 001 

r~ 

0. 001 

0. 0001 

0. 0001 

0. 0001 


lo"^ 

10-5 

10-5 

10"^ 

10-5 

10-5 

1 

o 

1. 1 X lo"^ 

1. 3 X 10"® 

4. 8 X lO"*^ 

1. 1 X lO"® 

1. 2 X lO"® 

4, 0 X lO"*^ 

,o-« 

1. 8 X lo"^ 

2.1x10"® 

3. 7 X 10"® 

l.|9 X lo”^ 

2. 1 X lO"® 

3. 6 X lO"® 

.0-5 

7. 9 X lO"'^ 

7.3 X 10"® 

1.0 xlO"® 

8.0 X lo"^ 

7.4x10"® 

1.0 X 10"® 

.0-5 

1.6 X 10“^ 

1.6 X 10"^ 

2. 0 x lO"® 

1.7 xlO"^ 

1.6x10"'* 

2.0x10"® 

.0-5 

2.6x10"^ 

3.4x10"'^ 

4. 0 X 10"^ 

2.6x10"^ 

3.4x10"^ 

4.0x10"® 


). Effect of Cq, D^, and t on the Position of Evaporation Boundary, y (cm) 



0.001 

0.001 

0.001 

0. 0001 

0.0001 

0. 0001 

; 

lO’^ 

10-5 

10-6 

10-" 

10-5 

10-6 

LO ^ 

5. 1 X lO"® 

5. 1 X lO"® 

■ i -8 

5.1 X 10 

6. 0 X lO"^ 

6. 0 X lO"^ 

6.0 X lO"^ 

10-’ 

5. 8 X lO"® 

5.19x10"® 

5.9xl0"® 

1.4x10"® 

1. 5 X lO'® 

1.4 X lO"® 

LO""^ : 

7.4x10"® 

7.4 X 10"® 

7. 5 X lo"® 

2.9 X lo"® 

2.9x10"® 

3*0x10'^® 

LO"'^ i 

8.4x10"® 

8. 5 X 10‘® 

8.7X10''® 

4,0 X l0“® 

4.0x10"® 

4.2xl0"® 

10-’* 

8.9x10"® 

9.0x10"® 

9. 2 X 10"® 

4. 4 X lO"® 

4. 5 X 10 ® 

4.8 x lO"® 
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It The position of the solidification boundary, z , increases with the evapora- 
tion-solidification time, t, but according to no apparently simple rules or , 
equations, at least from t = 0. 001 to 0. 02 second. Within this time interval, 
each increasing time step (which is 2. 0 or 2. 5 times the previous time t) 
increases the z value by multiplication factors F of 20. 1, 4. 3, 2. 2, and 
1. 6, respectively, from short to long t. (Fig. 6) 

2. Bj also exerts great influences on z (Fig. 6 and 7). Decreasing (due, 
for example, to low gravity) decreases z, by an average factor F of 8. 4 per 
decade of decrease in I>j . The average factor F is, however, greater for 

changing Dj from 10“^ to lO”^ cm^/sec (F = 8.8) than from lO”^ to 10 ^ 
cm2/sec (F = 7. 9) as shown in Fig. 6 and Table 10, 

3. The effect of on z depends very much on the evaporation-solidification 
time t. That is, and t strongly interact. Each decade increase in I>£ may 
increase z by as little as a niultiplication factor F of 1. 8 (at t = 0. 001 sec- 
ond and Bj from 10"® to 10"® cm^/sec), or as much as 10.6 (at t = 0. 005 
second and Bj from 10"® to 10"'^ cm2/sec), as shown in Table 10 and Fig. 6 


Table 10. Multiplication Factors F Due to Changing Bj at Different t 


t, sec 

0.001 

0. 002 

— 

0.005 

0.01 

0.02 

Average 

2 

Bj change, cm /sec 







lO"^ to lO"^ 

1.8 

5.4 

7.1 

7.8 

8.5 

7.9 

10 to 10 

6,4 

8.6 

10.6 

10.4 

7.7 

8.8 

Average 

4.1 

7.0 

8.9 

9.1 

8.1 

8.4 


4. There is a consistent trend of Cq on z also. This Cq effect on z is, how- 
ever, minimal. Each decade increase in Cq increases z by only about 3%. 
(Fig. 8) These z increases are practically zero at small Bj (10"® and 10“® 
cm^/sec) and low concentrations (Cq - 0. 0001). 

5. The solidification velocity z. changes with t in a manner stx’ongly dependent 
on D^ , not only in magnitude but in the shape of the z-t curve (Fig, 9). 
Such results are not available from close-form equations or simple com- 
puter simulations. 

Finding 4 indicates that it is impractical or Impossible to duplicate or study 
the effect of zero gravity on z, the position of the solidification boundary, by simply 
varjdng the antimony concentration (Cq) in germanium. If this trend holds generally, 
one may conclude that a single composition is probably adequate for an entire binary 
system, for the study of zero gravity effects on z. 
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Multiplification Factor in Z 


;^from 10 ' to 10 cm /sec 


“•6 "5 2 
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Fig. 6 Multiplication Factor for Crystal Length at Different Times for Two 
Decade^Chauges in ^ 
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Fig. 7 Position of Solidification Bo 
Mass Transfer Coefficients 
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Time, ms 

Fig. 8 Position of Solidification Boundary at Different Times for Various Initial 
Antimony Concentration In Ge 
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5.0 10.0 20.0 

ms 

Ltlon Velocity With D and t 


3 



Finding 3 suggests that an optimum time interval may exist (i.e. , from t * 0. 005 
to 0.01 second) for the study of effects on z of zero gravity through changing Dj , 
because of the large multiplication factors F in this time interval. Such study re- 
quires, however, submicron (depth) resolution techniques for sample characteriza- 
tion. 


For example, the effect of zero gravity is shown to reduce the effective liquid 
mass transfer from 10“^ to 10"5 cm^/sec (Ref. 10). Finding 3 thus suggests that 
?ero-gravity effects related to this mass transfer will be most noticeable when the 
crystal lengths on ground-base and space-processed samples are compared at 5 to 10 
ins after the initiation of cooling and solidification. Future space processing experi- 
ments may, therefore, be planned with this in mind, to maximize the chance of detect- 
ing and characterizing the generally subtle zero-gravity effects among the many in- 
herent noises. (Ref. 3). 

D. Summary of Combined Evaporation and Solidification 

The work under this contract has accomplished the following: 

1, Developed computer programs to study combined evaporation and solidifica- 
tion with surface evaporation and radiative cooling. 

2, Confirmed the importance of evaporation, as pointed out in the final report 
of our last NASA contract (Ref. 3), 

3, Using the developed computer programs, studied the individual effects of 
concentration on combined evaporation-solidification in the Sb-Ge system. 

4, Further studied the combined effects of effective liquid mass transfer and 
melt composition on combined evaporation-solidification through systematic 
computer simulation, 

5, Discovered the complicated interaction of melt composition Cq, liquid mass 
transfer coefficient Dy , and cooling time t on evaporation-solidification. 

6, Examined these new results in relation to possible zero-gravity effects and 
space processing. 

Notice that the computer programs developed \mder this contract, while still 
limited to a unidirectional model, accoxmt for not only heat and mass transfer, but 
also for surface evaporation and radiative cooling. Some predicted results have 
been confirmed, and some other unexpected results were discovered. For example, 
it is shown that the solidification velocity does not necessarily change monotonic^ly 
with the cooling time. In fact, the same velocity may increase or decrease with time, 
or reach maximum at an intermediate time, depending on D^ , Also, Dj markedly 
affects the position of the solidification boundary z, but not that of the evaporation 
boundary. On the other hand, increasing Cq affects the position of the evaporation 
boimdary, but not that of the solidification boundary (at least for Ge-Sb alloys studied). 
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The discovery of a maximum effect of a given change in , as measured by z, 
is probably the most directly relevant to space processing. By utilizing this informa- 
tion, many future space processes can be systematically designed for maximum sci- 
entific return. Such designs should involve not only the sample configuration and size 
or furnace capability, but also sample characterization techniques. The results on 
the position of the evaporation boundary, y (see Table 9) have also been analyzed. 
Among the three studied parameters, , Cq, and t, Cq most influences y. De- 
creasing Co by one decade from 0. 01 and 0. 001 reduces y by factors 8. 5 and 2. 6, 
respectively, as shown in Fig. 10. Increasing t always increases y, but at de- 
creasing rates with increasing t. This time effect is most noticeable at low Cq 
(F ig. 6). , on the other hand, is imimportant on y, particularly at high D/ values 

(Fig. 6). 

The evaporation velocity, y, decreases with increasing t and is reduced by 2-4 
orders of magnitude within 20 ms, depending on Cq. There is also an Interaction be- 
tween Cq and t on y. Thus, while a diffusion-controlled, somewhat steady-state 
condition is reached at low concentrations (c© = 0.001 and 0.0001) within 1 ms, the 
same condition is not reached at a higher concentration (Cq = 0. 01) after 2 ms (See 
Fig. 12). 

Thus, while greatly affects the solidification boimdary position, z (or single 
crystal length) or growth velocily z (Figs. 6, 7, 9), it has almost negligible effect on 
y (Fig. 11). This can be explained as follows: 'At reduced (such as in space), en- 
hanced solute pile-up ahead of the solid-liquid interface causes lowering of the melting 
point and, hence, reduced solidification under given thermal conditions. Because at 
970OC antimony evaporates seven to eight orders of magnitude faster than germaniiun, 
the solute depletion at the evaporating boundary is so rapid that even one or two orders 
of magnitude variations in change y only very slightly. 

E. Additional computed results 

Additional computed results for combined evaporation and solidification, and 
their uses in space processing, are given in Appendix E on an "Econotechnical Model 
of Capital Growth", 
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Evaporated Distance 






Time, ms 


Fig, 11 Position of Evaporation Boundary at Blfferent Times for Various 
Liquid Mass Transfer Coei^iclents 





4. THE GaAs STATISXrCALL Y DESIGNED EXPERIMENT 

The GaAs designed experiment was set up to achieve the following objectives:,^ 

1. Provide samples for checking our numerical computation programs. 

2. Study the effect of melt and growth variables on solute segregation and crystal 
quality, 

3. Show the power of such designed experiments, particularly in space processing 
where the sample numbers are necessarily limited and the cost of experiments 
is too high to tolerate inefficiency. 

A. Experimental Design 

Altogether 16 crystals were scheduled to be grown, according to a statistical de- 
sign. 

Four input or independent variables were studied. They are: 

1. Dopant types, at four levels: four in number: 

a. Cr 

b. Si 

c. Te 

d. Zn. 

2. Dopant concentrations at four levels: 

17 

a. 10 atoms/cc 

17 

b. 5 X 10 atoms/cc 

18 

c. 10 atoms/cc 

d. 5 X 10^^ atoms/cc. 

3. Growth rates, also at four levels: / 

a. 0. 10 in/hr 

1 i ■ . ../I':.-' 

b. 0. 16 in/hr 

c. 0. 22 in/hr 

d. 0.28 in/hr 
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4. Temperature gradient, also at four levels: 

a. 4° C/in. 

b. 6° C/in. 

c. 8*^/in. 

d. 10° C/in. 

Table 11 gives the GaAs crystal growth plan. 

The dependent variables or crystal characteristics to be measured at four dif- 
ferent positions along the crystal axis (i. e. , at 1, 1^, and 2 inches from the seeding 
planes) are: 

a. Free carrier concentration, in atoms/cc 

b. Hall mobility, in cm /volt-sec 

c. Resistivity, in ohm-cm. 

Zinc is a p-type dopant in GaAs, while Cr, Si, and Te are n-type dopants. Thus, 
the experiment gave both p and n-type GaAs crystals, with holes and electrons, re- 
spectively, as the majority carriers. Electrons have much higher mobilities than 
holes, as will be shown in this experiment also. 

j^The doping concentration in this experiment varies over a factor of 50, from 1 
X 10 to 5 X lOl® atoms/cc. This covers a range wide enough to study, e. g, , the 
effect of dopant concentration on segregation patterns or carrier mobility, without 
being too wide to cause difficulties in crystal growth or characterization. 

Growth rate and liquid temperature gradient are two well-known growth param- 
eters. Both affect not oiiy the growth morphology but also the solute distribution or 
segregation pattern. 

Among the dependent variables to be measured, the free carrier concentrations 
tie in with the solute redistribution pattern. Without contamination, compensation, but 
with 100% dopant ionization, these concentrations give the solute concentration profiles 
directly. The Hall mobility and resistivity are important characteristics for electronic J 
materials such as GaAs. 

' The experimental plan was designed according to a completely orthogonal Latin 
Square Design (Ref. 11). This design allows the study of a maximum of five input 
variables at four levels each. In this experiment, however, only four input variables 
were included. The design Is given in Table 11. 


4-2 



Table 11* GaAs Crystal Growth Schedule 


Crystal 

No. 

Dopant 

Concentration 

Growth 
Rate (in. /hr. ) 

Temp. 

Grad. (°C/in. ) 

1 

Te 

1 X 10^*^ 

0.16 

8 

2 

Si 

5 X 10^® 

0.22 

8 

3 

Cr 

5x10^® 

0. 28 

6 

4 

Si 

5 X 10^"^ 

0. 16 

6 

5 

Zn 

- , nlS 

5 X 10 

18 

0. 16 

4 

6 

Cr 

IxlOi 

0. 16 

10 

7 

Te 

lx 10^® 

0. 22 

6 

8 

Zn 

1 X 10^® 

0.28 

8 

9 

Cr 

5 X 10^"^ 

0. 10 

8 

10 

Si 

lx 10^® 

0.10 

4 

11 

Si 

1 x 10^*^ 

0. 28 ■ 

10 

12 

Zn 

1 X 10^^ 

0. 10 

6 

13 

Te 

5 X 10^*^ 

0. 28 

4 

14 

Zn 

SxlO^"^ 

0.22 

10 

15 

Cr 

1 X 10^"^ 

0.22 

4 

16 

Te 

5 x 10^® 

0.10 

10 


B. Experimental Detail 


The GaAs crystals were grown by Materials Research Corporation at Orangeburg, 
New York under contract. Furnace 3 was used exclusively for this experiment. This 
is a horizontal growth furnace measuring about five feet long and encased in a six-foot 
hood provided with access doors on all sides. A 2-inch (O. D. ) quartz tube about 2 
feet long is half filled with polycrystalline GaAs material together with the necessary 
dopants, and a single crystalline seed at one end. The growth furnace has three zones 
whose temperatures can be individually controlled to achieve the necessary temperature 
gradients in the liquid melt. After melting and homogenization and seed wetting, the 
furnace was moved away from the seed and growth began. Single crystals were not 
always obtained and some tests had to be repeated. Each crystal grown was later 
sliced by four diamond cuts at 1 , T, and 2 inches from the seed. 

During the crystal growth, it Was imexpectedly found that some vibrations en- 
hanced the chances of obtaining single crystal throughout the entire crystal length. 
These vibration effects have been systematically studied and will be discussed in the 
following section. 
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C. Experimental Data 


Macrophotographs on the three samples on each crystal were first taken. Ex- 
amples of such photos are shown in Fig. 13, for Crystal No. 1. This crystal is thus 
seen to vary in width from 33 mm at i inch, through 39. 5 mm at 1 inch and 42. 3 mm 
at li inches to 44. 3 mm at 2 inches. 

Data on free carrier concentrations. Hall mobilities, and resistivites at 1, 
and 2 inches from the seed are given in Tables 12, 13 and 14. 


Table 12. Free Carrier Concentration 



Intended 

Carrier 







Coef. 

Varia. 

Cone 

Avg/ 

intend 


Crys. 

No. 

Concentration 

(lOW) 

at 
2 in 

at 

1 in. 

at 

1| in. 

at 

2 in. 

Avg. 

Stand. 

Dev. 

m 

1 

1. 00 

2. 00 

2.00 

2. 09 

2. 13 

2. 06 

EH 


2. 06 

1.065 

2 


54.9 

53.1 

51.7 

50. 6 

52.6 



1.10 

0.922 

3 


44.5 

46.8 

47.7 

48.9 

47.0 

0. 21 


0.89 

1.099 

4 

50.0 

37.0 

39.0 

39.5 

39. 9 

38.8 

0. 14 

0. 036 

0. 74 

1.078 

5 

50.0 

48.8 

49.7 

49.9 

51.0 

49.8 

0.11 


0.97 

1. 045 

6 

10.0 

11.0 

11.0 

11.4 

12.7 

11.4 

0.11 


1.10 

1.155 

7 

10.0 

17.1 

20.3 

21.9 

23.8 

20.8 

0.32 


1.71 

1.392 

8 

10.0 

16.8 

17.7 

18.3 

20.4 

18.3 

0.17 

o.Oio 

1.68 

1. 214 

9 

5.00 

3. 84 

3. 97 

4. 13 

4.30 

4.06 

0.22 

0.055 

0.77 

1.120 

10 

10. 0 

15.6 

15.8 

16.7 

17. 8 

16.5 

0.11 

0.006 

1.56 

1. 141 

11 

12* 

1.00 

1.00 

1.47 

1.48 

1.62 

1.77 

1.58 

0. 15 

0.092 

1.47 

1.204 

13 

14* 

5.00 
5. 00 

6. 15 . 

6.60 

6.74 

6. 98 

6.62 

0.40 

0.061 

1. 23 

1. 135 

is 

1.00 

1.56 

1.72 

1.84 

1.89 

1. 75 

0.16 

0. 092 

1.56 

1.212 

16 

50.0 

57.1 

56.3 

56.1 



54.9 

56.1 

0.11 

0. 002 

1. 14 

0. 961 

♦Crystals 14 and 12 were not available. 


Avg. 


0. 030 

1.28 

1. 125 











at 1 1/2 Inches 


at 2 inches 


Fig. 13 Macrophotos on Sample Slices 1-1, 1-2, and 1-3 for GaAs Crystal 
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Table 13 Hall Mobility Data 


Crys. 

No. 



3* 

2 , 

^cm /v-sec) 

Stand. 

Coef. 

Ratio 


at 1 in. 


at 2 in. 

Avg. 

Dev. 

Var. 

2 in/i in. 

1 

2156 

2153 

2141 

2136 

2146. 5 

9.71 

0. 00453 

0. 991 

2 

1371 

1402 

1409 

1429 

1402. 8 

28. 17 

0. 02008 

1.042 

3 

237(i 

2362 

2357 

2342 

2357. 8 

13. 60 

0. 00577 

0.988 

4 

3330 

3327 

3320 

3318 

3323. 8 

5. 83 

0. 00175 

0.996 

5 

163 

164 

167 

160 

165. 8 

2. 43 

0. 01465 

1. 037 

6 

2376 

2365 

2366 


2366. 8 

7. 77 

0. 00328 

0. 993 

7 

3247 

3252 

3242 

3217 

3264. 5 

17. 00 

0. 00521 

0. 991 

8 

138 

144 

147 

149 

144.5 

5. 34 

0. 03697 

1.080 

9 

2279 

2264 

2250 

2239 

2258. 0 

19. 43 

0. 00860 

0. 982 

10 

2010 

2007 

2007 

1988 

2003. 0 

10.68 

0. 00533 

0. 990 

11 

12* 

3401 

3399 

3374 

3368 

3385. 5 

16.03 

0. 00473 

0. 990 

13 

14* 

2719 

2706 

2701 

2693 

2704. 8 

12.63 

0. 00467 

' 

0. 990 

15 

1480 

1476 

1465 

1462 

1470. 8 

8.74 

0. 00594 

0.988 

16 

2438 

2430 

2427 

2401 

2424. 0 

17.97 

0. 00242 

0. 985 


* Crystals 12 and 14 are not available. 
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Table 14 Resistivity Data (In 10 ohm-cm) 


Crys. 

— 





Stand. 

Coef. 

Ratio 

No. 

at j in, 

I' 

at 1 in. 

at l| in. 

at 2 in. 

Avg. 

Dev. 

Var. 

2 in. in. 

1 

14. 48 

14.47 

14.41 

14, 31 

14. 42 

0.083 

0. 0057 

0. 988 

2 

83.0 

94.0 

96.0 

97.0 

93.0 

7. 77 

0. 0835 

1.193 

3 

1.9 

1. 85 

1. 76 

1. 64 

1.79 

0.126 

0. 0705 

0. 863 

4 

5. 09 

5.16 

5.19 

5.27 

5.18 

0.087 

0.0169 

1.035 

5 

7. 85 

7. 63 

7.41 

t.41 

7.58 

0.214 

0.0282 

0.944 

6 

240. 

230. 

210. 

170. 

212. 5 

34.0 

0.160 

0.708 

7i 

1.09 

1.34 

1. 57 

k. 03 

1.51 

0.457 

0.302 

1.862 

8 

26.9 

25.2 

1 

23.7 

21.6 

j 

24.35 

2. 57 

0. 106 

0.803 

9 

7.15 

7.16 

7.10 

7.06 

7.12 

0. 0486 

0.0068 

0. 987 

10 

2. 00 

1.85 

1. 84 

1.74 

1.86 

0.126 

0.0168 

0. 870 

11 

12.5 

12.7 

12.9 

13.4 

12.88 

0.437 

0. 034 

1.072 

12* 









13 

3.06 

3.18 

3. 29 

3.36 

3.22 

0.146 

0. 0452 

1.098 

14* 









15 

270. 

260. 

220. 

190. 

235. 

38.9 

0.165 

0. 704 

16 

4.49 

4.32 

4.31 

4.33 

4.36 

0.087 

0. 0200 

0.964 

Avg 

1. 006 

*Crystals 12 and 14 were 

not available. 






D. Results and Conclusions 


Analysis of the experimental data (Ref. 11) have yielded the following results 
1. On Free Carrier Concentrations 
a. Concentration Variations Along the Crystal Length 

1, Without exceptions, this concentration changed monotonicaily along the crystal 
length. 

2. Concentrations increased in 12 of the 14 crystals. Only in No. 2 and 16 did 
the concentration decrease with distance from the seed. For the 12 crystals, 
the average ratio of concentration at 2 inches to that at ^ inch was 1. 176 (or 
17.|6% increase in concentration from ^ inch to 2 inches). The same average 
ratio was 0. 942 (or 5. 8% decrease in concentration from | inch to 2 inches) 
for the other two crystals. The combined average for this ratio on all 14 
crystals was 1.125. Thus, on the average, the concentration increased from 
I inch to 2 inches by about 12, 5%. These small increases indicate that 
nearly steady-state conditions were reached after I inch. 
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3. Cr-doped crystals had the greatest increases in concentration from J inch to 
2 inches by an average of 14. 6%; Te-doped crystals next, by 13. 8%; followed 
by Zn-doped crystals (12. 6%) and Si crystals (8. 6%). 

4. The maximum ratio of concentration at 2 inches to that at | inch (last column 
of Table 12) was 1. 392 (or 39. 2% increase from J inch to 2 inches) on Te- 
doped crystal No. 7; the minimum ratio was 0. 922 (or 7. 8% decrease from | 
inch to 2 inches) on Si-doped crystal No. 2. 

5. Increasing the growth rate seemed to increase this ratio, from 1. 074 at 
growth rate of 0. 10 in. /hr to 1. 163 at growth rate of 0. 28 in. Ar at least 
during the initial transient conditions. This conclusion agrees with the com- 
mon solidification theory based on the solute -enriched layer ahead of the 
solid-liquid interface (Ref 12). 

6. No systematic effects of melt concentration or liquid temperature gradient 
have been observed. 

b. Ratio of Average Concentration to Initial Melt Concentration 

These ratios are given in Column 10 of Table 12; 

1. The average concentrations in the crystals were generally higher than the 
concentrations in the initial melts, by about 28% on the average. 

2. More crystals (10) had higher average concentrations presumably due' to 
contamination; fewer crystals (4) had lower average concentrations prob- 
ably due to compensation, than the initial melt compositions. 

3. Maximum ratios of 2. 00 and 1. 71 (100% and 71% higher) occurred on Te-doped 
crystals No. 1 and 7; minimum ratio of 0. 74 (26% less) occurred on Si- 
doped crystal No. 4. 

4. On the average, Te-doped crystals had the highest ratios, averaging 1. 52, 
Cr-doped crystals the lowest with average 1. 08, Si (1. 22 average) and Zn 
(1. 33 average) crystals had intermediate ratios. 

5. Increasing the growth rate increased this ratio, from 1. 16 at 0. 10 in. /hr, 
through 1, 20 and 1. 46 at 0. 16 and 0. 22 in. /hr respectively, to 1. 32 at 0. 28 
in. /hr. 

C. On the Coefficient of Variation on Concentration 

This coefficient is also a measure of the variation of carrier concentration 
along the crystal length, due to the combined effect of doping, contamination» and 
compensation. 
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1. The average coefficient of variation on carrier concentration was 0. 030 


2. Again, doping types had important effect. Cr-doped crystals had the highest 
coefficients, averaging 0. 040. Si-doped crystals were next, averaging 

0. 034, followed by Te-doped crystals, averaging 0. 028. Zn-doped had the 
lowest coefficients, i. e. , the least variable, averaging 0. 006. 

3. Doping concentrations had the most important effect on these coefficients; 
increasing the initial melt concentration reduced the coefficients. The 
averages were: 0. 072 at lO^^ atoms/cc; 0. 051 at 5 x lO^*^ atoms/cc; 0. 010 
at 10^® atoms/cc, and 0.003 at 5 x 10^® atoms/cc. 

4. Increasing growth rate also increased these coefficients. Averages were: 

0. 020 at 0. 10 to 0. 16 in. /hr; 0. 37 at 0. 22 in. /hr, and 0. 42 at 0. 28 in. /hr. 

2. On Mobilities 

a. Average Mobility 


1. The average mobility of all 56 readings (four on each crystal) was 1587 
cm2/v-sec. 

2. Dopant type again had the most important effect. Zn-doping gave the lowest 
mobility. The average was only 155 cm^/v-sec, compared to the average of 
2340 cm^/v-sec for the other three types of doping. The ratio of 15. 1 
(2340/155) agrees with the fact that the Zn-doped crystals are p-type 
materials with holes as the majority carriers, while the other three types 
are all n-type materials with electrons as the majority carriers. Electrons 
move much faster than holes in semiconductors such as GaAs. 

3. Among crystals other than Zn-doped, Te-doping gave the highest mobility 
(2604 cm^/v-sec) and Cr-doping the lowest (2075 cm2 /v- sec), with Si- 
doping in the middle (2371 cm^/v-sec). 

4. As expected, increasing the doping concentration in the initial growth melt re- 

duced the average mobility, from 2465 cm2/v-sec at initial melt concentration 
Coof lO^*^ or 5^ lol*^ atoms/cc, through 1223 cm^/v-sec at c = 10^® atoms/ 
cc, to 1074 cm2/v-sec at 5 X lOl® atoms/cc. ° 


b. Mobility Variation Along Crystal Length 


1. With only one notable exception on Crystal No. 7, mobility changed monotoni- 
cally along the crystal len^h. 


2. Eleven crystals decreased, and only three crystals (No. 2, 6, and 8) 

Increased in mobility from the seed ends on. This is most easily seen by 
examining the ratio of mobility at 2 inches to that at i inch. The average ratio 
for the 11 crystals was 0. 989 (or 1. 1 % decrease from | inch to 2 inches). The 
average ratio for the 3 crystals was 1. 053 (or 5. 3% increase from inch to 
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2 inches). The combined average for all 14 crystals was 1. 003 (or only 

0. 3% increase from | inch to 2 inches). 

3. The maximum ratio was 1. 080 (or _8. 0 % increase in mobility from | inch to 
2 inches) on Crystal No. 8, and the minimum ratio was 0. 982 (or 1. 8% de- 
crease in mobility from | inch to 2 inches) on Crystal No. 9. 

4. Increasing the growth rate v increased this ratio, from 0. 985 at v = 0. 10 
in. /hr, through 1. 004 at v = 0. 16 in. /hr and 1. 007 at v = 0. 22 in. /hr, to 
1. 012 at V = 0. 28 in. /hr. 

5. Increasing initial doping concentration in the melt also increased this ratio, 
averages being 0. 990 at c = 10^^ or 5 x 10^"^ atoms/cc to 1. 013 at c - 10l8 
or 5 X 1018 atoms/cc. 

c. On Coefficient of Variation on Mobility 

This coefficient is also a measure of the variation of mobility along the crystal 
length, due to the combined effect of doping and crystal perfection. 

1. The average coefficient of variation on mobility for all 14 crystals was 0. 0062 
(or 0. 62%). Maximum coefficient was 3. 70% with Crystal No. 8, and the 
minimum coefficient was 0. 18% with Crystal No. 4. 

2. Zn-doping gave decidedly more mobility variation along the crystal length 
than the other three types of doping. Averages were 0. 0233 (or 2. 33%) vs 

0. 0050 (or 0. 50%). 

3. Increasing the growth rate v increased the mobility variation along the crystal 
length. Averages were 0. 0043 (or 0. 43%) for v = 0. 10 and 0. 16 in. /hr vs 

0. 0084 (or 0. 84%) for v = 0. 22 and 0. 28 in. /hr. 

4. Increasing doping concentration in the initial melt also increased the mobil- 
ity variations aloi^ the crystal length. Averages were 0. 0046 (or 0.46%) for 
Cq = 10^*^ and 5 x 10^^ atoms/cc vs 0. 0083 (or 0. 83%) for Cq - 10^® and 

5 X 10^® atoms/cc. 

3. On Average Resistivities 

1. The mininuun resistivity of all 56 readings (four on each crystal) was 
1. 64 X 10"®, and the maximum 0. 27 ohm-cm, 

2. Among the four different types of dopants, Cr gave the highest average 
resistivity (28. 2 x 10"® ohm-cm), while Te the lowest (4.2 x 10"® ohm-cm), 
with Si (10. 4 X 10”3 ohm-Cm) and Zn (12. 3 x 10”3 ohm-cm) giving intermedi- 
ate values. This conclusion ^rees with the common observation that Cr is a 
’*semi-insultaing" dopant in the GaAs while Te can give very low resistivity 
materials. 

3. As expected, the higher the dopant concentration, the lower the resistivity. 

Averages are 35.2 x 10"3 4,91 x 10“3, 11,0 x 10”3, 8.6 x 10-3 ohm-cm 

at concentrations Cq = 10^ \ 5 x lO^'^, 10^®, and 5 x 10^® atoms/cc. 
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E. Discussion on the GaAs Experiment 

Crystal growth is generally still an art. Also, two of the 16 crystals originally 
planned in this experiment could not be grown because of unexpected personnel changes 
at our v^iidor. The measurement errors appeared to be large, when one compares 
the measured results given in Tables 12 to 14 with other measured results independently 
obtained by Prof. Franklin Wang’s group at the State University of New York at 
Stony Brook (See Table 15). 

Yet the results appear to be useM. In particular, the conclusions seem to check 
well with established practices in the industry (Ref. 13). There are no anomalies 
or unique results observed. 

The ratios of concentration at 2 inches to that at 5 inch give some idea about the 
segregation coefficients. The higher this ratio, the greater the segregation co- 
efficient. Since Te-rdoped crystals had the highest ratios, 1, 52 on the average, hence, 
Te would liave the highest segregation coefficient in the GaAs melt and Cr the lowest 
segregation coefficient in the GaAs melt. Si and Zn would have intermediate 
segregation coefficients, with the segregation coefficient of Zn being greater than that 
of Si. Exact segregation coefficients on all these four dopants in GaAs are not avail- 
able in the literature. 

One interesting confirmation of the existing theory (Ref. 12) is the conclusion 
that increasing the growth rate increased the above ratio. This can be explained as 
follows; increasing the growth rate increases the solute pile-up ahead of the solid- 
liquid Interface, at least for the beginning portions of completed growth. The localized 
liquid solute concentration at the interface is increased. The solute concentration 
in the solid crystal freezing out also is increased to maintain a nearly constant segrega- 
tion coefficient for these dilute melts, hence, the observed results. 

One disturbing thing is that the melts in some cases were contaminated or 
compensated. Some crystals had the solute concentrations decreasing along the 
length Instead of increasing. The sources of this contamination is not known. 

Contamination of semiconductor melts is not uncommon (Ref. 13) in view of the 
extreme dilution of the doping elements. Generally, the solute concentrations of the 
solutes varied in this experiment from 10^*^ to 5 x 10^® atoms/cc, corresponding to 
ppm or tens of ppm. A single dust particle near the seed woidd thus easily have 
contaminated the melt and even compensated the tra<?e doping elements. 

The contamination problems, however, make the strict correlation with 
computed profiles difficult. There are also only four data points on each crystal even 
though we tried hard but failed to develop other more or less Continuous probing 
techniques to obtain more data. 

Nevertheless, we believe that this experiment fulfilled two of the stated objec- 
tives of the experiment, i. e. , to study the effect of melt and growth variables on 
solute segregation and crystal quality; and to show the usefulness of such designed 
experiments. No correlation between computed and measured results has been made, 
partly because of budgetary constraints. 
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In particular, we have shown that this type of statistically designed experiment 
works, in crystal growth also. The good correlation of the observed results with the 
generally accepted results lend credence to this approach. We also feel that this type 
of experiment As particularly useful in space. This is because space experiments 
must necessarily be limited in sample size, and yet zero-gravity effects may be so 
subtle relative to the many ever-existing noises (experimental errors or other effects) 
(Hef. 3) that improved techniques in experimental design and data analyses must be 
resorted to for meaningful results to be obtained with a minimum number of samples, 
in a short timet and at low costs. 


Table 15 Data on GaAs Samples, Work Done at SUNY 


Crys. 

Resistivity, ohm-cm 

2 * 

Mobility, cm /volt, sec 

Concentra- 
tion, lO^"^ 

1 

0.00193 

2212 

14. 6 

2 

0. 00529 

981 

12.0 

3 

0. 00565 

2160 

5. 12 

3 

0.00533 

2296 

/ 

5. 11 

4 

0. 0026 

1122 

21.4 

4 

0.00232 

1649 

16.3 

5 

0.0348 

90 

19.9 

5 

0. 032 

84 

23.0 

6 

0.170 

2082 

0.17 

6 

0.0756 

( 

5325 

0. 155 1 

7 

0.00192 ' 

1232 

26.4 ! 

1 

7 

0.00192 

1265 

25.7 

8 

0. 0357 

90 

19.4 

8 

0.0339 

120 

15. 3 

10 

0.00335 

1611 

11.6 

10 

0,00287 

1455 

14. 9 

11 

0.0166 

2840 

1.32 

11 

0.0156 

2722 

1.47 

13 

0.018 

2463 

1.40 

13 

0. 0176 

2987 

1. 19 

15 

0. 594 

989 

0.106 

15 

0.755 

1021 

0,0811 

16 

0.00208 

1213 

24.8 
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5. SEGREGATION IN GaAs CRYSTALS 


In connection with our statistically designed GaAs experiment, we selected a 
number of GaAs (k in) slabs for microprobe analyses. Studies on the Ga and As con- 
centration profiles reveal an unanticipated segregation phenomenon and some evidence 
of strong convection currents. This section details the experimental procedures, re- 
sults, data analyses, conclusions, and discussions. 

A. Experimental 

The GaAs crystals were grown at our vendor, Materials Research Corporation 
at Orangeburg, New York. The growth equipment was horizontal Bridgeman quartz 
tube furnace. The GaAs polycrystalline melt materials and a single crystalline 
seed were placed in a sealed 2 inches diam x 2 feet quartz tube or boat filled with 
argon imder partial vacuum. Only about half of the sealed quartz tube was filled with 
GaAs. 


The sealed quartz tube was placed inside the quartz tube, which was supported 
in a SiC tube about ^ inch- thick. The SiC tube is heated by three sets of equally 
spaced Globars inside an insulated furnace (Fig. 14). The furnace is thus longi- 
tudinally divided into three individually controlled temperature zones. Specifically, 
each zone can be controlled in temperature to about independently of the other 
two zones. The furnace is driven by a scan motor to move longitudinally relative to 
the SiC and quartz tubes containing the sealed quartz tube and GaAs melt for the 
growth to initiate and continue. 

It was found that various kinds of vibrations existed on the growth furnace. 
Further, the vibrations from the scan motor markedly and unejq)ectedly improved the 
crystal quality and yield. Details of these vibration effects will be described in Sec- 
tion 6. 

17 18 

The GaAs were usually doped with 10 to 5 x 10 atoms/cc of Cr, Si, Te, and 
Zn. Typical growth conditions were 0. 19 in. /hr with a liquid temperature gradient 
of70C/in. 

After growth, the crystals were sliced on an ID diamond saw into ^ inch slices. 
Selected slice samples were mounted and microprobed at 3-micron intervals up to 
72 microns with an electron microprobe Model AMR/3 made by Hiillips. 


B. Results 

On each GaAs slice sample, a number of microprobe transverses were taken. 
Figure 15 shows locations of the surface points from which concentration profiles 
perpendicular to the surface were made on a typical sample. 

The Ga and As concentration profiles related to these surface points were given 
in Fig. 16 to 21. 
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Fig. 14 GaAs Growth Furnace 
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Fig. 15 Location of the Ga and As Microprobe Transverses 
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Fig. 17 Ga and As Concentration Profiles at Point 2 
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Fig. 19 Ga and As Concentration Profiles at Point 4 
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C. Analysis of Results 


At all the six surface locations, Ga is depleted, while As is enriched. The 
distance, d, over which this depletion or enrichment takes place varies from 15 
microns at location 1, to 47 microns at location 4, The values of d are given in Table 
16, line 2 . 

Beyond this distance d, the Ga and As conce7i|?rations in the bulk of the crystal 
are constant, down to the third decimal place; Ga at 48, 3 weight percent (50. 1 atomic 
percent), and As at 51. 7 weight percent (49.9 atomic percent). Notice that the sum 
of Ga and As concentrations add up to 1, 000, as they should. The same is true for any 
probed point on the sample. 

The surface concentrations of Ga and As, however, vary from location to loca- 
tion. Ga concentration varies from 0,411 at surface point 2 to 0,458 at surface point 
6, (line 7 of Table 16) while As concentration varies from 0. 542 at surface point 6 to 

0. 589 at surface point 2 (line 4 of Table 16). Averse surface concentrations are Ga 
0.4422 and As 0. 5578 weight fractions (or 46.0 and 54,0 atomic percent). These 
surface concentrations are also given in Table 16. 

Table 16 (lines 6 and 9) also gives the average concentrations (of Ga and As) for 
each transverse, and the ”h^ distance" d/2 showing the depths inside the surafce 
that have these average concentrations. 

D. Discussion 

The precision of the microprobing is better than 0. 001 in weight fractions, for 
both Ga and As. This can be seen from the following: 

1. Underneath the surface regions (i.e, , 15 to 47 microns thick), the Ga con- 
centration is constant at 0.483 ± 0.000; and As concentration is constant at 
0. 517 ± 0. 000, for the (up to 18 times at location 1) repeated readings at 3 
micron intervals. 

2. The sum of Ga and As concentrations at any point on the sample adds up to 

1. 000 ± 0 , 000 . 

The observed concentration variations must, therefore, be real. 

These surface concentrations cannot be caused directly by surface evaporations, 
because these evaporations are limited generally to the fractional micron surface lay- 
ers, as shown in Section 2. 

Also, if the segregations were due to evaporation, we would except the surface 
regions to be depleted in As and enriched in Ga, because of vapor pressure consider- 
ations given below. Ihis is opposite to what was observed. Table 17, taken from 
p. 57 of Ref. 14, for example, gives the temperatures (in ®C) for As and Ga to reach 
different vapor pressures (in atmospheres). Notice the much hi^er Ga evaporating 
temperature than As for the same vapor pressure. Also, the density, molecular 
weight, and mole fraction of Ga in the evaporating GaAs are roughly equal to the cor- 
responding values for As (Ref. 13). 

We must, then, assume that the observed Ga and As segregation patterns are 
caused by solidification. 
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Table 16. Characterization of GaAs Sample ^ 


Loca- 

tion 

1 

2 

3 

4 

5 

6 

Average 

Line 

1 

d 

15 

Distances 
45 45 

47 

38 

19 

34.8 microns 

2 

d/2 

4.0 

8.9 

6. 5 

7.0 

6.9 

5.4 

6.45 microns 

3 



As concentrations (Weight Fractions) 



Surface 

0.551 

0. 589 

0.546 

0.545 

0. 574 

0. 542 

0. 5578 

4 

Interior 

0. 517 

0. 517 

0. 517 

0. 517 

0. 517 

0.517 

0. 5170 

5 

Average 0. 534 

0.553 

0. 532 

0. 531 

0. 546 

0> 530 

0. 5377 

6 



Ga Concentrations (Weight Fractions) 



Surface 

0.449 

0.411 

0.454 

0.455 

0.426 

0.458 

0.4422 

7 

Interior 

0.483 

0.483 

0.483 

0.483 

0.483 

0.483 

0.4830 

8 

Average 0.466 

0. 447 

0.468 

0.469 

0.454 

0.470 

0.4623 

9 



Effective Segregation Coefficients (Surface/interior) 


As 

1.066 

1.139 

1.056 

1.054 

1.110 

1. 048 

1. 079 

10 

Ga 

0.930 

0.851 

0. 940 

0.942 

0.882 

0.948 

0.9155 

11 




Estimated Growth Velocities 




0.775 

0.180 

0.537 

0.530 

0.278 

0„ 739 

0.5065 cm/sec 12 


Table 17. Evaporating Temperatures To Reach Given Vapor Pressures 


Pressures 

0. 0001 

0.001 

0.01 

0.1 

0.5 

1.0 

Atm 

As 

308 

363 

428 

499 

578 

610 

°C 

Ga 

1178 

1329 

1515 

1751 

1985 

2071 
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The exact leveling off of the Ga and As concentrations beyond the surface regions 
indicates the attainment of steady-state conditions. The surface segregations are, 
therefore, merely transient solidification results. 

The surface regions are very thin (15 to 47 microns) relative to the size of the 
samples (oblong shape 24 mm x 46 mm). Hence, the initial melt concentrations also 
were Ga 0. 483 and As 0. 517. 

The effective segregation coefficients for Ga and As in GaAs melt are, therefore, 
about 0.9155 and 1.079, respectively. However, these coefficients, obtained from the 
surface and initial concentrations, vary from location to location, as given in Table 16 
(lines 10 and 11), 

Notice that the effective segregation coefficient for Ga in GaAs varies from 
0. 851 at location 2 to 0. 948 at location 6, a range of 11. 4%. The same coefficient for 
As in GaAs varies from 1. 048 at location 6 to 1. 139 at location 2, a range of 8. 0%. 

The surface concentrations, Cf, of As (or Ga) appeared to depend on the half 
distance, d/2 as shown in Fig. 22. A regression analysis (Ref. 15) shows, however, 
that the correlation is statistically insignificant. For example, the regression coef- 
ficient is only 0. 6918, the value of t is 1.916, and the value of F is only 3.672. All 
these values are below the 90% significance levels (Ref. 15). Thus, the surface con- 
centration Cj is independent of the half distance d/2, 

The surface region thickness dg also appears to depend on the half distance d/2, 
as shown in Fig. 23, A regression analysis confirms this. The (regression) coef- 
ficient is R = 0. 8426, t = 3. 129, and F = 9. 791. These are all significant at the 95% 
levels (Ref, 15). 

The prediction equation is: 

tg = -11.95 + 7.253 X d/2 ±8. 55, t^ and d/2 both in microns, 
where ±8. 55 is the adjusted standard error of estimate (Refs. 16-17). 

From the surface and bvilk (interior) concentrations, and the half distances, one 
can estimate the local freezing velocities V, following the Pohl diffusion-limited 
model (Ref, 18). These V values are given in Table 16, line 12. 

As indicated above, the microprobing accuracy is good, probably to better than 
0. 0005 weight fraction. The observed Ga and As concentration variations are, there- 
fore, real. 

Therefore we can safely say that Ga was depleted from the surface and As en- 
riched at the surface. The primary mechanisms of this depletion and enrichment 
could not be evaporation. They must be the result of solidification. 

That is, Ga segregated from the GaAs melt having an initial melt concentration 
of 48. 3% by weight, with an effective segregation coefficient less than unity. The 
average coefficient was 0. 9155. On the other hand. As segregated from the same 
GaAs melt having an initial melt concentration of 0, 517, with an effective segregation 
coefficient greater than unity. The average coefficient was 1. 079. 




Fig. 22 Relation Between Maximum As Concentrations and Half Distance, 
or Effective Segregation Coefficients 
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Fig. 23 Correlation Between Surface Region Thiclmess and Half Distance 
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According to Hansen’s phase diagram (Fig. 24) for GaAs (Ref, 19), the liquidus 
curve is a single continuous line extending across the entire composition range. The 
tangent to the liquidus curve at the GaAs compound composition is a nearly horizontal 
line. There is no solidus curves except for the vertical line at the GaAs composition 
(See Fig. 24). There should, therefore, be no segregation due to solidification. The 
solid and liquid phases must be of the same composition, i.e. , 50 atomic percent of 
Ga and 50 atomic percent of As. (or 0. 482 weight fraction of Ga and 0. 518 weight 
fraction of As). 

On the basis of our segregation results, the actual phase diagram for GaAs at the 
peak region may look like that shown in Fig. 25. The initial GaAs melt may be slightly 
poorer in Ga than 50 atomic percent, e.g. , 0.481 weight fraction, not as observed. 
Possibly, the probing accuracy is off by a substantial calibration error, even though 
the precision was good. The GaAs freezing out will thus follow the solidus line PG, 
and not PA. Ga will then be surface depleted, as observed. But the bottom of the 
solidus line PG must then be near 50 atomic percent (0. 483 weight fraction of) Ga. 

The pealc point P must thus be slightly off to the right of the 50:50 (atomic percent) 
mark. 


Also, both the effective segregation coefficients (for Ga and As) are close to 
unity. Hence, the liquidus and solidus lines for either side of the GaAs peak must be 
very close at the peak, since the liquidus at the peak is nearly horizontal, the solidus 
lines must also be nearly horizontal. There must then be a wide spread between the 
two solidus lines in Fig. 25. This is also not present in the Hansen diagram (Fig. 24). 

The observed As concentration was as high as 58. 9 weight percent (or 57. 2 
atomic percent). The deviation of the solidus line PG at close to the melting point of 
GaAs (1511®K according to Ref. 20) must be at least 7. 2 atomic percent, or almost 
one major grid spacing on the Hansen's phase diagram (Fig. 24). Such is also not 
shown. 

We conclude from the above that Hansen’ s GaAs phase diagram needs modifica- 
tion, to account for the observed solidification phenoroen. 

Let us examine more closely the relationships between the surface concentration 
surface region thickness dg» half distance d/2, and effective segregation coeffi- 
cients kg. The regression between cf and d/2 is statistically insignificant. The re- 
gression of d and d/2 is significant at the 95 % level. The regression of kg for Ga or 
As on Cf is significant at the 99% level. 

Thus, while C£ and d/2 or d can vary independently of each other, d and d/2 can 
be estimated one from the other. At any location for either Ga or As, kg can be di- 
rectly computed with high confidence, once the surface concentration cf is known by 
probing. 

The growth velocities v from Pohl analysis (Ref. 18) show a unique trend. In 
particular, location 2 had the smallest velocity, followed lay that at location 5, while 
locations 1 and 6 had the highest velocities. The high velocities at locations 1 and 6 
can be explained on the basis of cooling effects of the contacting bottom at location 6, 
and the corner effect at location 1. But the unique velocity distribution for the other 
locations, 2 to 5, still deserves careful interpretation. 
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Let us first review briefly the assumptions involved in the Pohl analysis (Ref. 

18). Pohl assumes that the one-dimensional solidification is diffusion-limited, and 
that the segregation coefficient is not veiy large. The effect of limited stirring (due 
to convection, for example) can in some sense be lumped into an overall effective 
liquid mass transfer, whose practical results follow the Pick’s first law (Ref. 21). 

The effective segregation coefficient, ke, in this case is, however, close to unity, 
quite different from the assumed values. But the Pohl equation is still, in this partic- 
ular case, better than the Tiller equation (Ref. 12) which is good only for ke approach- 
ing zero. 

In spite of these limitations, we feel that the Pohl analysis still should give some 
rough estimates of the localized growth velocities. The unique velocity distribution 
pattern must, therefore, be significant. 

This unique velocity pattern can not be explained directly on the basis of thermal 
phenomena, such as temperature fluctuations, because then the velocities on the top 
surface locations should at the worst change monotonically. Such was not the case. 

The temperature inside the sealed quartz boat or tube was in fact very stable, because 
of the enclosing quartz and SiC tubes, and equidistant, spaced-apart Globar heating 
rods outside the SiC tubes. . The thermal inertia of the SiC and quartz tubes and the 
surroimding ambient should have amply stabilized the temperature Inside the sealed 
quartz boat or tube. The three- zone, temperature controllers also controlled the 
temperature of each zone to within about 0. 5®C of their intended values (Ref. 22). 

The unique velocity distribution could not have been caused by physical or chemi- 
cal processes related to the melt surface, such as oxidation, chemical reactions, or 
dust settling, because the velocities should then be uniform or vary monotonically 
from the center location 5 to the edge or corner location 1. The effect of extraneous 
vibrations observed (and to be described in Section 6) also could not be the cause, for 
similar reasons. 

By this process of elimination, we thus conclude that the unique velocity distri- 
bution must be related to phenomena inside the GaAs melt, such as convection currents, 
as shown below. 

When the furnace was slowly withdrawn, the SiC tube cools first, followed by 
the quartz tube, the sealed quartz tube or boat, and the GaAs melt. With temperature 
controls on the SiC tube the cooling must be very uniform. 

The bottom (location 5) of the GaAs melt, which contacted the sealed boat, must 
cool fairly fast, as was observed. The corner (location 1) might have cooled equally 
fast, or even slightly faster, because of the combined effects of heat radiation and 
comer cooling. The locations 5 and 2 were probably subject to high normal fluxes of 
convection currents, and received heat from the interior of the hot GaAs melt. These 
two locations therefore cooled the slowest. This leaves the intermediate locations 4 
and 3 to cool at intermediate speeds, as observed. 

If this explanation proves to be correct, its implication to containerless, space 
processing is clear. Specifically, in zero-gravity conditions, there are no gravita- 
tional convection currents; only true diffusion operates. For a spherical geometry, 
there will only be a radial temperature gradient and concentration gradient. The 
growth velocities will be the same at all surface locations (for a simple isometry). 
There will be no variations in temperature, concentration, and growth velocity in 
directions other than radial. This definitely was not the case with our GaAs crystal 
growth conditions, as can be seen by comparing Fig. 16 to 21. 
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Not only will be the space GaAs crystals be more imiform in concentration, but 
the crj-^stals will be more perfect. This is because the crystal can tolerate uniform 
radial shrinking stresses, but not nohradial stresses. There are several reasons for 
this : (1) the neighboring regions stressing each other were at different temperatures 
and, therefore, had different strengths; (2) absence of arch effects for nonradial 
stresses, and (3) higher temperature and concentration gradients introduces higher 
mismatch stresses. 

Note that GaAs has great potential as an electronic material, but is generally 
found to be of lower quality than expected. The mobilities, for example, are very 
low, particularly in p-type materials. Usually, the mobilities in these p-tjq>e mate- 
rials are only between 100 to 500 cm2/volt-sec (Ref. 13). These poor character- 
istics are generally believed to be caused by imperfections. We do not know how 
much of these imperfections are related to convection currents. But certainly, re- 
moval of the convection-induced defects must help. 

In other words, space processing, because of the absence of convection cur- 
rents, mi^t give better GaAs crystals in terms of uniformity and defect contact. 

E. Summary of Conclusions 

This section describes some results of microprobing work on GaAs crystal sam- 
ples. We find that; 

1. Probing gave Ga and As concentrations precise to ±0.0005 in weight fractions. 

2. The bulk of the original GaAs melt had a Ga concentration of 0.483 weight 
fraction (50. 09 atomic percent) and As concentration of 0. 517 weight frac- 
tion (49,91 atomic percent). 

3. There was significant Ga depletion and As enrichment on the surfaces of from 
15 to 47 microns thick. Average surface concentrations were: Ga - 44.22 
wei^t percent 54.00 atomic percent) and As - 55.78 wei^t percent 

(46. 00 atomic %). 

4. These segregations are real, because of the small probing errors, 

5. These segregations could not be caused by evaporation, but must be caused 
by solidification. 

6. On solidification, Ga segregated with an average effective segregation coef- 
ficient of 0. 9155, and As with an average effective segregation coefficients of 
1. 079. 

7. These segregation coefficients kg are not anticipated, at least according to 
Hansen's phase diagram on GaAs, 

8. The half distances, or depth to average concentrations, varied from 4, 0 to 
8, 9 microns, with an average of 6.45 microns. 

9. The surface concentration (of either Ga or As) did not correlate with the half 
distances d/2. 
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10. The surface region thickness dg correlated significantly (95%) with, and could 
be predicted from, d/2. The prediction equation is: 

dg (microns) = -11. 95 + 7. 253 x d/2 (microns) ±8. 55 

11. The effective segregation coefficients can almost be exactly calculated 
from the surface concentrations alone: 

ke for As = 1. 936 

kg for Ga = -0. 084 - 2. 064 

12. The growth velocities at different locations have been estimated by a PoW. 
analysis. These velocities varied from 0.180 to 0. 775 cm/sec. These initial 
velocities apply only to the surface regions (15 to 47 microns thick). 

13. The melt interface probably was concave toward the liquid. 

14. Hansen's GaAs phase diagram probably needs modification. 

15. These growth velocities have a urique distribution pattern which can not be 
explained on the basis of furnace temperature, surface reactions, dust 
settling, vibration, or other causes external to the GaAs melt. The pres- 
ence of a gravitational convection current, however, quite satisfactorily ex- 
plains this velocity distribution pattern. 

16. In space, there is little gravitational convection current. There should be 
no temperature, concentration and stress gradients in the nonradial direc- 
tions. The stress pattern will be relatively simple and readily tolerated. 
More uniform and defect-free GaAs crystals will, therefore, be grown in 
space. 
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6. VIBRATION EFFECTS ON CRYSTAL GROWTH 


Our vendor, the Materials Research Corporation of Orangeburg, New York, in 
growing the 14 GaAs crystals for our statistically designed experiment, noted that a 
damaged and unbalanced scan motor (used for moving the horizontal Bridgeman Furnace) actu 
ally improved the yield of single crystals. "Yield" as used here means the proportion 
in a growth run of crystals which are completely single crystalline over the entire 
length of about 2 feet. The location of this motor on the furnace is shown in Fig. 26. 


That the yield was improved was convincingly shown by distinct improvements 
whenever this motor was switched back and forth from one furnace to another of iden- 
tical design. Data comparing the yield with the unbalanced motor to the yield with a 
balanced motor as control are given in Table 18. These data on GaAs crystal growth 
(Ref. 22) are for two identivally designed furnaces during three successive growth 
periods, or runs 1 to 3, after the motor was unbalanced due to accidental damage and 
the unusual effects of this damaged motor were noted. For example, during growth 
period 1, furnace 3 was with the unbalanced motor and had three completely single 
crystals out of four grown. Simultaneously, furnace 10 had the balanced motor and 
yielded no completely single crystals out of six tried. 

Table 18. GaAs Yield Data Comparing Motors, Furnaces, and Growth Periods 


Growth 

Period 

Furnace 3 

Furnace 10 

1 

Unbalanced motor; yield 3/4 

Balanced motor; yield 0/6 

2 

Balanced motor; yield 0/3 

Unbalanced motor; yield 3/3 

3 

Unbalanced motor; yield 3/3 

Balanced motor; yield 1/5 


It can be shown that such results are statistically significant (Ref. 23). That is 
the unbalanced motor markedly Improved the yield of completely single- crystalline 
GaAs crystals, on both furnace 3 and furnace 10, for all periods 1 to 3. 


Systematic work on the effect of vibration during crystal growth is seldom docu- 
mented and studied in detail. According to Dr. R. A. Laudise, Assistant Director, 
Materials Research Laboratory of the Bell Telephone Laboratories, no systematic 
work has been done on the effects of vibration on crystal growth, nor does he know of 
any current work in this area at other laboratories (Ref. 24). hi some previous MIT 
work (Ref. 25),, vibrations of varying frequencies were introduced to the crystal grower, 
without any appreciable effect being noticed at any frequency. 

A. Relation of Vibration Effects to Low-G Ciy stal Growth 

Vibrations may have beneficial or harmful effects on crystal yield. We hypothe- 
size that the vibration effects we observed probably have to do with liquid mixing or 
mass transfer ahead of the solid-liquid interface. Some Skylab experiments also 
pointed to reduced liquid mixing and effective solute transfer as a most important 
zero-gravity segregation effect, according to E. McKannon, A. Ukanwa, and 
T. Bannister (Ref. 10). We have, therefore, taken a close look at vibration effects 
during GaAs crystal growth. 
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We arrived at the above hypothesis through a process of elimination, even though de- 
tailed thermal, metallurgical, mechanical, and thermodynamic analyses have not been 
made. 


A notable possible effect of vibration to crystal growth is to change the instan- 
taneous growth rate by, e. g. , momentarily perturbing the temperature or temper- 
ature gradient through adiabatic expansion and cooling at the solid-liquid interface. 
Another possible effect of vibration is to modify the instantaneous solute concentra- 
tion through localized and periodically varying (non- symmetric), artificial gravity 
effects due to density differences between the solute and solvent melt. When the in- 
stantaneous growth velocity, temperature gradient, or solute concentration, sepa- 
rately or in combination, exceeds the critical limits, constitutional supercooling and 
polycrystalline growth result. But since the vibrations from the unbalanced motor is 
only of the milli-g level, as will be shown, they should have relatively negligible effects, 
particularly when compared to the effects of unavoidable, fractional degree temper- 
ature fluctuations always present in the melt. There may also be a boimdaiy layer 
that damp these effects through mechanical or thermal inertia effects. For similar 
reasons, the periodic tensile and compressive or shear forces at the interface can be 
also ruled out. It has been estimated that the effect of applying 100, 000 psi in some 
respects is equal to a change in melt temperature of only 0. 3®C (Refs. 26-28). Such 
temperature fluctuations are inherent in the furnaces. 

Vibration may also force the liquid melt to make better contact with the wall of 
the quartz tube, at least for ultrasonic welding or soldering operations. But the fre- 
quency and the milli-g forces are again probably too low to be very effective. In addi- 
tion, this effect of vibration will introduce additional nuclei and is opposite to what 
happened. Another possible vibration effect is to precipitate nuclei in the GaAs melts 
saturated or nearly saturated with some unknown impurities. But again this should re- 
sult in nucleation, and not denucleation, as was probably observed. 

We can also assume that the GaAs melt contained refractory, particulates that 
could act as nuclei for polycrystailine growth. Such particles may, for example, be 
dust particles or grains from the quartz tube (Ref 29). These particles must, in 
addition, be difficult to dissolve in, or wet by, the GaAs melt. Vibrations, however, 
enhanced the wetting and dissolution, thereby achieving denucleation (Ref. 30). How- 
ever, the milli-g forces at the particular frequencies are probably too small or incor- 
rect to achieve these results. 

The most plausible mechanism, we believe, is as follows: Vibration introduced 
some additional mixing effects in the liquid ahead of the interface. A controlled, 
nlilli-g, vibrative force might have introduced a gentle "massaging" action into the 
liquid. The longitudinal •‘Vibrations (along the axis of the quartz growth tube) increased 
the longitudinal, effective liquid mass transfer and reduced the likelihood of longi- 
tudinal, constitutional supercooling. Transverse vibrations (perpendicular to the 
growth tube axis) increased the transverse, effective liquid mass transfer and reduced 
transverse solute concentration gradients and constitutional supercoolings. In either 
case, polycrystalline growth is supressed, as was observed. 

Crystal growth is generally still an art. Growth defects often appear and dis- 
appear mysteriously. In particular, many crystals imcontrollably go polycrystalline. 
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Fig. 26 Schematic of Typical Furnace 
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Fig. 27 Instrumentation for Eecording Furnace Vibrations 
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Yet without good crystals, modem electronics would not have been developed and 
many other physical and metallurgical discoveries would not have been made. Many 
electrical, optical, and mechanical structures are still material limited. B would 
thus appear that efforts should be spent to follow up on significant leads related to 
crystal growth and quality improvement. In addition, the space environment may 
provide iinique conditions that may lead to improved growth techniques or theoretical 
understandings, as will be shown in the discussion ot' this vibration work. We have 
undertaken a study to determine exactly what happened due to the vibrations from the 
unbalanced motor. 

B. Experimental Detail 

We simultaneously measured the vibrations in the horizontal and vertical direc- 
tions with Glennite AUR 250-30G accelerometers, together with Model 8300 Crescent 
amplifiers and an oscilloscope for instant visual observations. The amplified vibra- 
tion signals were recorded with a Honeywell Model 5600 tape recorder on l/2 inch 
magnetic tapes operating at 60 inches/sec ui the FM mode. Figure 27 shows the in- 
strument set-up. 

Five critical points in the furnace were selected for the vertical and horizontal 
vibration measurements. These points are given below and are marked "X" in Fig. 28. 

1. On scan motor 

2. On aluminum chain block supporting the lead screw 

3. On SiC tube supporting the quartz tube 

4. hi quartz boat inside the quarlz tube at the seed location 

5. On heating oven for quartz tube. 

The Glennite accelerometers together with the Crescent amplifiers have a com- 
bined sensitivity of about 1. 5 volts DC per g (gravitational constant), and a linearity of 
db 0. 25% of full scale. The frequency response of tiie amplifiers is -3 dB at 300 Hz and 
-12. 5 dB at 103 Hz. The noise content in the output is 10 mV which is equivalent to 
0.0069 g. 

Under the operating conditions given above, the bandwidth of the recorder is lO'^ 
Hz with a signal to noise ratio of 49 dB. The noise output of the recorder is equivalent 
to 0. 0035 g. An additional channel in the recorder was used for voice identification of 
each run which insured positive identification of the data during playback. An oscil- 
loscope, used to monitor the signals being recorded, identified any overloads or 
malfunctions of the equipment durii^ recording, and provided a check on the overall 
input-output characteristics of the recorder. 

Before making a tape recording, either furnace 3 or furnace 10 and the recording 
point on it were selected. In all runs furnace 3 had the unbalanced scan motor, while 
furnace 10 had the balanced or control motor. The unbalanced motor could simply be 
turned on or off. The balanced motor, on the other hand, could not only be turned on 
or off, but also be adjusted to different voltage settings to change its rotating speed. 
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Furnace 3 had the cooling fan for the electronic temperature circuit in operation, but 
furnace 10 had the same cooling fan removed for repair during the time of measure- 
ment. 

The background noises were first recorded. The several motors on the furnace 
or in the room were then (selectively or) successively turned on. In general, the 
order of turning on the motors were; cooling fan, heater fan (for room heating), hood 
fan (for furnace exhaust), and scan motor (for furnace driving relative to the stationary 
quartz boat containing the GaAs melt for growth). Each motor was turned on for 30 
seconds before the next motor was activated, except for the scan motors, which were 
turned on for one minute intervals. After the last or scan motor was on fer one minute, 
it was turned off, followed by the successive turning off of the other motors in the re- 
verse order at 30 second intervals. 

In Run 6, the end doors of the furnace hood (Fig. 28) were tapped with the middle 
finger, successively but in order, at locations marked "X” in Fig. 28. Some time 
was also spent waiting for a freight train to pass (within 200 feet). This allowed us to 
record the vibrations due to the passing train in Run 11. 

C. Data and Data Analyses 

1. Data; 

Altogether, 16 recording runs were made on four reels of magnetic tape. The 
contents of the recorded tapes are given in Table 19. 

2. Oscilloscope Pictures; 

The oscilloscope photos were made by playing back the magnetic tapes on the 
Honeywell Model 5600 tape recorder connected to an oscilloscope (Tektronix 545). The 
yertical and horizontal signals were Individually amplified before being fed siniulta- 
heously to the oscilloscope. Each oscilloscope photo thus contains two signals in the 
form of two wavy lines. The top line is for the vertical signal, while the lower one 
for the horizontal signal. On those photographs having 1/9 g per cm display, the noise 
component due to the instrumentation is about 0. 07 cm, or 0.0008 g. 

Figures 29a through 29e, from Run 7, are oscilloscope photos of the vibration 
signals inside the quartz tube (at the seeding plane) of furnace 3. They are, respec- 
tively, for backgroimd noise, cooling fan on, heater fan on, hood fan on, and scan 
motor on. The hood fan introduced some vibrations, mostly in the vertical direction 
(0.06 g). The unbalanced scan motor, however, introduced both vertical and horizontal 
vibrations of far greater magnitude (0.10 g). 

Figures 30a through 30d, from Run 8, are oscilloscope photos of the vibration 
signals inside the quartz tube but on furnace 10. They are, respectively, for back- 
ground noise, heater fan on, hood fan on, and (balanced) scan motor on. The balanced 
scan motor did not seem to introduce any significant vibrations, in either the vertical 
or horizontal direction. 

Figures 31a through 31d, from Run 5» are oscilloscope photos of tlie vibration 
signals on the SiG tube on furnace 3, respectively for background with cooling fan on, 
heater fan on, hood fan on, and scan motor on. By comparison with Figures 29a-29e, 



Table 19 Vibration Signals on GaAs Furnaces as Recorded 


Reel Run Furnace Record 
No. No. No. Point 


1 1 

1 2 

1 3 

1 4 

1 5 


3 Sic tube 
3 SiC tube 
3 SiG tube 
3 SiC tube 
3 SiC tube 


2 6 


3 


SiC tube 


2 7 


3 


Quartz boat 


2 


8 


10 


Quartz boat 


3 9 3 Sic tube 


3 10 3 Oil Scan motor 


Notes 

Background- -30 secs recorded 

Heater fan only on — 30 sec 

Hood fan only on — 30 sec 

Scan motor and hood fan on — 30 sec 

Cooling fan on — 30 sec 

Heater fan on — 30 sec 

Hood fan on — 30 sec 

Scan motor on — 60 sec 

Scaii motor off— 30 sec 

Hood fan off — 30 sec 

Heater fan off — 30 sec 

Cooling fan off — 30 sec 

Tapping three hood doors. 

Side doors 9 times each, front door 20 times 

Background— 30 sec 

Cooling fan on — 30 sec 

Heater fan on — 30 sec 

Hood fan on — 30 sec 

Scan motor on — 60 sec 

Scan motor off — 30 sec 

Hood fan off — 30 sec 

Heater fan off— 30 sec 

Cooling fan off — 30 sec 

Background — 30 sec 

Heater fan on — 30 sec 

Hood fan on — 60 sec 

Scan motor at 5 mv — 30 sec 

Scan motor at 10 mv — 30 sec 

Scan motor at 15 mv — 30 sec 

Scan motor off — 30 sec 

Hood fan off — 30 sec 

Heater fan off — 30 sec 

Background — 30 sec 
Cooling fan on — 30 sec 
Heater fan on — 30 sec 
Hood fan on — 30 sec 
Scan motor on— 60 sec 
Scan motor off — 30 sec 
Hood fan off— 30 sec 
Heater fan off— 30 sec 
Cooling fan off — 30 sec 
Background— 30 sec 
Cooling fan on — 30 sec 
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Table 19 Vibration Signals on GaAs Furnaces as Recorded (Continued) 

Reel Run Furnace Record 

No. No. No. Point Notes 

Scan motor on — 60 sec 
Scan motor off— 30 sec 
Cooling fan off — 30 sec 
Passing freight train within 200 feet 
Background — 30 sec 
Cooling fan on — 30 sec 
Scan motor on — 60 sec 
Scan motor off — 30 sec 
Cooling fan off — 30 sec 

Background — 30 sec 
Cooling fan on — 30 sec 
Scan motor on — 30 sec 
Scan motor off— 60 sec 
Cooling fen off — 30 sec 
Background — 30 sec 
Cooling fan on — 30 sec 
Scan motor on — 60 sec 
Scan motor off — 30 sec 
Cooling fan off — 30 sec 
Background— 30 sec 
Cooling fen on — 30 sec 
Scan motor on— 60 sec 
Scan motor off — 30 sec 
Cooling fen off— 30 sec 

Figures 33a and 33b, from Runs 12 and 10, show respectively the vibrations on 
the aluminum, support chain block and directly on the scan motor of furnace 3. By com- 
paring these two photos with Figs. 29e and 31d, one notices that, ^ expected, the vi- 
brations grew stronger as the recording point was shifted from inside the quartz tube 
(Fig. 29e), through the SiC tube (Fig. 31d) and chain block (Fig. 33a), to directly on 
the scan motor (Fig. 33b). Tapped end doors vibrated even more strongly (Fig. 29f). 

Figure 34, from Run 11, shows that no additional noises were introduced on 
chain block of furnace 3 by a freight train passing within 200 feet. 

D. Computer Signal Analyses 

The recorded vibration signals on the magnetic tapes were played back onto an 
Adage Ambilog computer for Fourier analyses of the spectrum distribution, power 
spectral density (PSD), autocorrelation, and cross-correlation. The Honeywell tape 
recorder was again used, to play back the signals at a rate of 3 x 10^ samples per 
second for storage in a buffer storage unit of 4 x 10^ capacity. When the data buffers 
were filled, each channel was unbiased by subtracting its means and then unitized by 
dividing by its maximum value. The PSD function of each series was then calculated 
and stored in the spectrum buffer, and subsequently read out and automatically plotted 


3 11 3 Chain block 

3 12 3 Chain block 


4 13 3 Chain block 


4 14 3 Sic tube 


4 15 3 On furnace 
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Fig. 29 V'ibrations on SiC Tube of Furnance No. 3 


•3-9 












6-10 







Fig. 31 Vibrations on SiC Tube of Furnance 3 
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i 


f 


the horizontal background noise appeared to be slightly increased. But the greatest 
differences between these sets of pictures are observed in the vibrations Introduced 
by the unbalanced scan motor. The amplitude in the horizontal direction is almost 
doubled. The vertical vibration signal also increased in amplitudes, particularly in 
the secondary vibrations of higher frequencies. 

Figures 32a through 32c, from Run 6, show that tapping the end doors at loca- 
tions marked as "X" in Fig. 28 introduced severe vibrations on the SiC tube on fur- 
nace 3. Note that the vertical scale on these photos is 0. Ill g/cm, rather than 
0.037 g/cm as in all the other pictures (Figs. 29, 30, 31 and 32). Also, the horizontal 
vibration signals were about 2-1/2 times stronger than the vertical signals. Figure 
32a also show's some transient behavior of the vibration in the horizontal direction. 



30 X gain 
5V/cm 
1.5 u/G 


a - Tap End Door 
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b - Tap End Door c - Tap End Door 
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Fig. 32 Vibrations Due to Tapoing on Sic Tulx? of Furnace 3 
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a - Vibrations on Chain Block 
Fig. 33 Vibrations Due to Scan Motor of Furnace 3 



Fig. 34 Vibrations Due to Passing Freight Train on Chain Block of Furnace 3 
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by a digital plotter. Figures 35 to 56 are examples of the PSD analyses of some of 
the vibrations recorded. The vertical ordinate in all these figures is proportional to 
the squared amplitude of the signal (e, g. , g2) while the horizontal abscissa indicates 
the frequency of the signal. 

By putting both square and sine waves of three different amplitudes each, we 
proved that the vertical ordinate on the PSD curves are proportional to the applied 
signals (or voltages) squared. Since we know the applied voltage is proportional to 
the acceleration force , (in g's), the vertical ordinate can thus be calibrated into g^. 
This has been done for each PSD curve, as is shown in the figures. 

Figures 35-44 show, in subgroups of two for the vertical and horizontal signals, 
the computer-plotted, PSD analyses of the vibration signals recorded on the quartz 
boat inside furnace 3. These figures are, in order: backgroimd noise; cooling fan on; 
heater and cooling fans on; hood, heater, and cooling fans on; and scan motor, and 
hood, heater, and cooling fans on. The background noises peaked at about 30 Hz in 
the vertical directions (Fig. 35), while the horizontal signal peaked at about 55 Hz (Fig. 
36). The cooling fan did not appreciably change the vibration spectra, either in the 
vertical or in the horizontal directions (Figs. 37 and 38). The heater fan did, however, 
introduce at least one new vibration spectrum with a peak at 50 Hz in the vertical direc- 
tion (Fig. 39), which is higher than the 30 Hz vertical background noise, but no appre- 
ciable horizontal peak (Fig. 40). Tlie hood fan noises peaked at 123 Hz, both in the ver- 
tical and horizontal directions (Figs. 41 and 42). The scan motor introduced at least 
two additional peaks, a strong peak at 54 Hz, which is much stronger than the hood or 
other noises, and a weaker peak at 247 Bjz (Figs. 43 and 44). 

Figures 45 to 53 show similar results of PSD analyses of the recorded vibration 
signals, but for furnace 10 with a balanced motor. The backgroimd noise here peaked 
at 120 Hz, both in the vertical and horizontal directions (Figs. 45 and 46). Neither 
the cooling fan (Fig. 47) nor the heater fan (Figs, 48 and 49) contributed much vibra- 
tion signal. The hood fan did, however, introduce at least one additional strong signal 
at 123 Hz,, with the vertical signal being much stronger than the background (Fig. 50), 
but the horizontal signal being weaker than the background (Fig. 51). The scan motor 
noises peaked at 28 Hz in both directions. (Fig. 52 and 53), 

Figure 54 shows the PSD analyses of the vibration signals on the SiC tube of fur- 
nace 3, with the scan motor, hood fan, heater fan, and cooling fan on. The signals of 
the scan motor again predominated at, as expected, 54 Hz. 

Figures 55 and 56 show the effect of tapping the end doors. The strong tapping 
signals peaked at 110 Hz, 

Table 20 summarizes the results of PSD analyses of the above vibration signals 
and gives such statistical data as average signal density, majcimum signal density, and 
standard deviation of the spectra for each case. The peak accelerations (in g) and the 
frequence of the peaks are also given. 

Figure 57 exemplifies an autocorrelation analysis of the vibration signal in the 
vertical direction, as recorded on the SiC tube of furnace 3, with the scan motor, hood 
fan, heater fan, and cooling fan on; 2, 000 data points were included in the analysis. 
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Fig. 37 PSD Inside Furnace 3, Vertical Position, Cooling Fan On 
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Source 


Position 


Furnace 


Direction 


X 


Bacl^ound 

Background 

Cool Fan 

Cool Fan 

Htr & Cool Fan 

Htr & Cool Fan 

Hd, Htr, Cool Fan 

Hd, Htr, Cool Fan 

t’ans, all 

Background 

Background 

Cool Fan 

Cool Fan 

Htr & Cool Fan 

Htr & Cool Fan 

Hd, Htr, Cool Fan 

Hd, Htr, Cool Fan 

Scan, Hd, Htr, Cool 

Fans, all 

Tapping On Door 
Tapping On Door 


Position 

Furnace 

Direction 

X 

X 

max 

a 

a 

max,g 

fa, 

*max 

Quartz bt 

3 

Vert 

-0. 272 

BB 

0.208 

0.016 

30 Hz 

<^artz bt 

3 

Hor 

-0.020 


0.294 

0.012 

120 

Quartz bt 

3 

Vert 

-0. 269 


0. 258 

0.016 

30 

Quartz bt 

3 

Hor 

-0.022 

-0.131 

0.135 

0.012 

120 

Quartz bt 

3 

Vert 

-0.272 

0.112 

0.196 

0.020 

130 

Quartz bt 

3 

Hor 

-0. 294 

0. 081 

0.265 

0.022 

125 

Quartz bt 

3 

Vert 

- .266 

-0.104 

0.237 

0.020 

125 

Quartz bt 

3 

Hor 

- .024 

- .088 

0.281 

0.020 

125 

Quartz bt 

3 

Vert 

-0. 272 

0. 139 

0. 272 

0.038 

57 

Quartz bt 

3 

Hor 

-0.027 

0.184 

0.356 

0.085 

57 

Quartz bt 

10 

Vert 

0.002 

-0.078 

0.307 

0.020 

25 

Quartz bt 

10 

Hor 

-0.078 

-0. 085 

0.201 

0. 012 

no 

Quartz bt 

10 

Vert 






Quartz bt 

10 

Hor 

0.510 


0.195 

0.010 

125 

Quartz bt 

10 

Vert 

0.002 


0.258 

0.020 

30 

Quartz bt 

10 

Hor 

0. 504 


0.218 

0. on 

no 

Quartz bt 

10 

Vert 

0.000 

-0.086 

0.284 

0.021 

35 

Quartz bt 

10 

Hor 

0.505 

0.078 

0.219 

0.011 

125 

Quartz bt 

10 

Vert 

0.009 

-0.075 

0.325 

0. 020 

28 

Quartz bt 

10 

Hor 

0.578 

0.077 

0.254 

0.012 

120 

SiC tube 

3 

Vert 
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0.220 

0.318 

0.022 

60 

SiC tube 

o 

_> 

Hor 

-0.167 

0.281 

0. 549 

0. 235 

60 

Sic tube 

3 

Vert 

-0.026 

0.193 

0.208 

0.018 

no 

Sic tube 

3 

Hor 

-0.130 

-0.328 

0.226 

0.042 

no 















Each of the 2000 points on this figure was obtained by correlating the signal spectra 
with a similar signal spectra but shifted by increments of one data point. The wavy 
curve fluctuates from +1.0 for perfect positive correlation to -1. 0 for perfect negative 
correlation. The regularity and periodicity of the wave line are to be noted, 

E. Other Analyses 


Other computer analyses show that the vibration spectra all appeared to be 
Gaussian type. Also, walldng and jumping near the furnace, or opening and slamming 
a door to the next room did not introduce any additional vibration. 

As expected, maximum vibration occurred directly on the scan motor, followed 
by positions on the supix)rt chain block (for the transmission gears), SiC tubes inside 
the furnace, and the quartz boat inside the furnace at the seed position. 

The imbalanced scan motor on furnace 3 generated substantial vibrations at 54 
Hz in both the vertical and horizontal directions. Finger tapping of the end doors also 
created transient vibrations whose frequencies between 60 and 100 Hz had the largest 
amplitudes in the horizontal direction. In the vertical direction, 80 to 120 Hz were 
the most predominant, hi both cases, 100 Hz signals had the largest amplitudes. By 
considering that tapping could excite all frequencies, including resonantly frequencies, 
these tests point towards a resonant frequency of about 100 Hz for the SiC tube in 
furnace 3. In amplitudes, these excited vibrations are twice that created by the im- 
balanced scan motor on furnace 3. 

F. Conclusions and Discussions on Vibration Effects on Crystal Growth 


We conclude, from the above results, that 

1. the fluctuations in crystal yield at our vendor were not due to random causes, 
but due to a definite causative effect. 

2. Such a causative effect can be traced to a damaged, unbalanced motor. The 
effect of this motor on crystal growth was significant, independent of the 
growth furnace or growth periods 

3. This motor introduced a idb ration spectra in the vertical and horizontal 
directions which were completely different from those due to an undamaged 
or control motor. 

4. The vibration spectra due to tins unbalanced motor has been fully character- 
ized by detailed measurements. Their amplitudes peaked at milli-g's at 

54 Hz. These amplitudes were much greater than the background noises. 

5. In contrast, the cooling fan for the electronic temperature control circuit did 
not Introduce any appreciable additional vibrations over the background; the 
heater fan for the room had some additional vibration spectra of low ampli- 
tudes at 50 Hz in the vertical direction. The hood exhaust fan introduced 
some vibrations at 123 Hz which were stronger than those due to the heater 
fan, but still much wealcer than those due to the unbalanced scan motor. 
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6. Tapping the end doors of the furnaces created vibrations at 110 Hz which v/ere 
much stronger than those due to the unbalanced motor. This indicates the 
need to improve the door designs to minimize vibrations, 

7. All the vibration spectra appeared to be Gaussian type. The regularity and 
periodicity of the wave line on the autocorrelation curve suggest the regu- 
larity and periodicity of the vibration signals. That is, these signals were 
probably sinusoidal in shape. 

8. Walking and jumping near the furnace, or opening and slamming a door (6 
feet away) to the next room did not introduce any additional measurable vi- 
brations. This conclusion confirms the good isolation of the furnace sup- 
port. 

9. A freight train passing within 200 feet did not noticeably affect any of the 
backgroimd vibration spectra. This disproves the common fallacy that good 
crystal growth facilities must be located* often at great costs, far away 
from railway tracks 

We also conclude that a proper vibration (or oscillating) force of even less than 
0. 05 g is significant in improving ciystal quality and yield. Such a small g force is 
subtle and is easily masked in ground experiments, such as the previous MIT work 
(Ref. 25), but probably requires space environment for systematic study, i.e., isola- 
tion, evaluation, and optimization. 

This vibration experiment presented some interesting questions related to space 
processing. Did the vibrations due to the unbalanced motor indeed improve the liquid 
mixing ahead of the solid- liquid interface? Did these ^dbrations really have a "mas- 
saging” action? If milli-g forces were beneficial, could proper micro-g forces also 
be beneficial? What are the vibration spectra due to g- jitters, astronauts walking or 
handling the equipment? Could some of the mysterious Sl^lab results be due to some 
hidden vibrations? 

We have thus carefully and systematically documented the different sources of 
vibrations on one set of GaAs growth furnaces, and, in particular, fiiUy character- 
ized the beneficial vibration spectra due to an unbalanced motor. We cannot, without 
any additional work, pinpoint exactly which particular single or combined vibration 
frequencies are beneficial, nor the critical amplitude ranges for each such frequency. 
Nor can we speculate as to the exact mechanism of improved crystal yield. This 
important crystal growth knowledge must probably also be unambiguously obtained by 
controlled space experiments for subsequent uses on earth. 
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7. CONCLUSIONS 


We have briefly reviewed results developed under our previous contract 
(NAS8-27891) and assessed where we were when present contract started. 

The new computer programs developed xmder this conti act are found to be 
powerful in the study of the combined effects of evaporation and solidification. Some 
unexpected results have thus been discovered. Specifically, we have found that the 
solidification-evaporation is sensitively affected in a complicated manner by the com- 
bined effects of melt temperature and composition, or liquid mass transfer mecha- 
nisms such as convection currents. The solidification velocity does not even neces- 
sarily change monotonically with the cooling time. Also, the liquid mass transfer 
markedly affects the position of the solidification boundary, but not that of the evapo- 
ration boundary. On the other hand, varying the initial solute concentration in the 
melt affects the position of the evaporation boundary, but not that of the solidification 
boundary. We thus have confirmed the importance of evaporation in solidification 
studies, and also showed some limitations of the conventional solidification theories 
(without evaporation). 

Also from these computed results, we found a maximum effect of a given change 
in the liquid mass transfer on solidification. By using this information, we may maxi- 
mize the subtle zero-gravity effects among the effects of other uncontrolled variables 
so that future space processing experiments can be designed systematically for maxi- 
mum scientific returns, particularly when the technique related to a proposed ''econo- 
technical Model of Crystal Growth" (Appendix E) is used. 

We also discovered, by microprobing, some imexpected segregation effects in 
GaAs samples, and revealed, through determining the velocity distribution pattern, 
some evidence of convection currents during the GaAs growth. We thus conclude that 
space-grown GaAs samples by the horizontal Bridgeman method should be more uni- 
form and perfect, in line with other Skylab crystal growth results. 

A GaAs statistically designed experiment is presented to show the efficiency 
and applicability of this type of experiments in crystal growth, particularly during 
space processing where the sample sizes and available equipment, personnel, and 
time are necessarily limited. 

In growing the GaAs crystals, we also noticed some unique, beneficial effects 
of vibration on crystal yield and perfection. We examined the many possible effects 
of vibration on crystal growth and believe that the most plausible mechanism relates 
to liquid mixing ahead of the solid-liquid interface. Hence, understanding vibration 
effects may contribute importantly to understanding the reduced liquid mass transfer 
in space conditions. We have carefully documented the different sources of vibrations 
on some GaAs growth furnaces and, in particular, fully characterized the beneficial 
vibration spectra that improved the crystal quality. However, more work is 
still required to confirm the proposed beneficial mechanism, and to optimize and 
apply these mechanisms in space and groimd crystal growth operations. 
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8. RECOMMENDATION FOR FUTURE WORK 


We recommend: 

1. Additional study of convection currents in crystal growth. 

2. Investigate the mechanism of the beneficial effects of vibration. 

3. Use of computed results to correlate and explain existing space processing 
results and to coordinate space crystal growth experiments, 

4. Continuation on the improvement of computer programs to develop additional 
crystal growth models. 

5. Additional systemmatic study of combined evaporation and solidification with 
improved programs. 

6. Study why segregation in GaAs and possibly other III-V compounds do not 
follow the usual phase diagrams. 
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APPENDIX A 


IMPORTANCE OF EVAPORATION 


Evaporation is important in space melting and solidification for the following 
reasons: 

• Significant evaporation of alloy components always occurs at high tempera- 
tures in space vacuum environments 

• High- temperature evaporation of alloys is generally a neglected area of 
systematic research. Yet, unless the complete evaporative segregation 
behavior is understood and analyzed, solidification and its related segregation 
effects may not be properly studied because of ill-defined initial conditions. 
Before the liquid alloy can be controllably solidified or even melted, there is 
invariably some surface evaporation to cause changes in composition, freezing 
temperature, supercooling characteristics, nucleation and growth morphology 
conations, and the like 

• Controlled space evaporation probably most closely meets the requirements 
of our model of normal evaporation. We may be able to obtain material 
purity or evaporation standards, thermal properties, or even such basic 
thermodynamic properties as heat of evaporation, activity coefficients, and 
sticking coefficients that are difficult or impossible to obtain on earth 

• Evaporation is a much simpler process than freezing, since the former does 
not involve such complicated phenomena as nucleation, phase transformation, 
and constitutional or nonconstitutional supercooling. Thus, in normal evapo- 
ration for specific geometries or alloy systems, we may ideally isolate and 
Investigate such other phenomena as heat conduction or radiation, liquid or 
solid diffusion, fluid dynamics, and convection currents. Exact knowledge of 
these phenomena is necessary to understand solidification 

• Evaporation causes surface cooling due to the heat of evaporation.. This 
evaporative cooling effect is particularly important in low-melting materials 
(Ref. Al) 

• Different rates of evaporation at various surface regions give rise to un- 
balanced forces and momenta that may produce erratic or imwanted accelera- 
tions, surface distortions and vibrations, exceedin^y large "equivalent 
gravities, " and possibly new types of powerful convection currents in zero- 
gravity conditions 

• Evaporation may cause the surface composition of certain unwanted or un- 
suspected impurities to be increased a thousandfold or millionfold within 
seconds so that the layer’s melting point and other thermophysical properties, 
nucleation characteristics, base for imdercooling, and critical velocity to 
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avoid constitutional supercooling may be completely unexpected. In fact, 
anomalous "constitutional" or evaporative melting on cooling, or solidifica- 
tion on heating, is possible because of surface evaporation. In addition, 
very large artificial gravities (e. g. , 10 g), strong fluid disturbances, or 
even new and significant convection currents may be produced from surface 
evaporation. These phenomena have been observed in the M553 movies, ac- 
cording to Dr. Martin Tobin of Westinghouse Co. , Pa, 

The much greater evaporative segregation effects, if unaccoimted for, would 
almost certainly conceal any minor or subtle zero-gravity effects, particularly in the 
presence of other unknown or uncontrolled effects. Definitive space solidification 
work should probably, therefore, be preceded by an evaporative compatibility study 
of the sample materials and their possible associated impurities. In fact, evaporation 
is almost certain to be very important or so overwhelming that the effect of zero- 
gravity or freezing segregation may be masked or even reversed. A freely sus- 
pended molten drop in space may, for example, have its surface solute concentration 
greatly enriched (as much as a millionfold), by neglected and undetectable trace im- 
putities within seconds of its deployment. We are then dealing at the critical surface 
layer with a completely new and unanticipated alloy having an entirely different com- 
position, melting point, surface tension, thermophysical properties, latent heat of 
fusion, undercooling and nucleation characteristics, growth morphology, and the like. 

From this we can also see that any anal 3 rtical, computer, or experimental study 
on solidification may yield completely unexpected or irrelevant results if the impor- 
tant and ever-present evaporation phenomena is not adequately taken care of. This is 
particularly true in the study of nucleation, undercooling, and space processing. 
Another important aspect of the present contractual work is to incorporate this gener- 
ally neglected evaporation phenomena to define the exact initial and boundary condi- 
tions before and during the alloy solidification process. 

References: 

Al. Reichman, J. "Solidification of Metal Spheres in Vacuum," Grumman Research 
Memorandum RM-544, June 1972. 



APPENDIX B 


COMPUTER PROGRAM WITHOUT SURFACE EVAPORATION AND RADIATION 


Solidification, even in one-g, is a complicated process involving a multitude of 
interrelated phenomena such as mass and heat transfer, phase change, and fluid 
motion. Comprehensive reviews on solidification have been given, for example, by 
Chalmers (Ref. Bl), Tiller (Ref. B2), Christian (Ref. B3), and Li (Ref. B4). 

Solidification in zero-g is still very complicated. Here, gravitational force is 
negligibly small, but other effects as a result become important. For example, sur- 
face tension often plays a dominant role in determining the sample shape, processing 
technique and the resulting contamination level of the processed samples. Evaporation 
is another ever-present, complicating or dominating factor, but one that may be used 
to advantage when understood. Neglected, or improperly controlled e^ auoration 
may drastically change sample surface composition, fluid motion, equivalent gravities, 
nucleation, and undercooling characteristics as previously described. The previous 
program, under Contract NAS 8-27891, however, does not deal with evaporation. 

Mathematical Definition of Solidification Problem 

To understand thoroughly solute segregation either from combined evaporation 
and solidification, or in single-crystal growth, one requires a complete characteriza- 
tion of the (mass) diffusion and temperature fields in the solid crystal and remaining 
melt. The zero-gravity effect on the solidification may be overshadowed by other 
effects invariably present (such as evaporation) in any such growth process — a con- 
dition necessitating that such characterization he accurately defined. Unfortunately, 
the coupled partial differential equations for the diffusion and tt mperature fields are 
generally not solvable. Although special case solutions have been given for some 
types of usually physically nonsatisfying, two-phase Stefan problems, for the general 
case solution we must resort to numerical computations. Existing numerical methods 
are always subject to such unrealistic assumptions as constancy of interfacial velocity, 
temperature or temperature gradients, segregation coefficients, diffusion constants, 
and other material thermo physical properties. 

Under NAS 8-27891, a number of computer programs were developed to study 
the unidirectional solidification of a binary alloy. These programs employ anal 3 d;ical 
and numerical methods. The analytic program is based on some closed-form solu- 
tions of a simple model and gives results for our numerical program to compare. The 
model for the analytic program deals with a binary alloy at a constant temperature and 
concentration throughout the initial liquid melt, with the surface temperature 
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instantaneously dropped below the liquidus temperature. The liquid-solid interface 
temperature is assumed constant, and the concentrations of the alloy at th? interface 
are given by the phase diagram having curved liquidus and solidus lines. In addition, 
the interface boundary plane moves according to a square root law relative to the 
solidification time. The program also allows the interface temperature and interface 
boundary to vary from these fixed rules, but in practice the variation is negligible 
and not above the computer error level (Ref. B5). 

Although covered in detail in the final report on NAS 8-27891, the mathematical 
formulation of the model is presented below for the sake of completeness. 


We deal in unidirectional solidification with a liquid binary alloy to be direction- 
ally solidified into two phases, liquid and solid. We consider the liquid alloy to be 
semi-infinite with original (at t = 0) temperature Tq and concentration C^. Solidifi- 
cation occurs when the temperature at x = 0 is changed from Tq to a lower value Tn , 
either instantaneously or gradually, so that Ti is below the temperature T 2 at whi^ 
the liquid mixture at concentration Cq can be in equilibrium with a solid phase. As 
solidification occurs, the solid phase grows and its boundary is located at x = y(t), and 
the interface temperature at this point is Ti(t). The partial differential equations de- 
scribing the solidification process are the following; 
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where the variables T, C represent the temperature and concentration (of solute in 
solvent) and the subscripts >£, s denote the liquid and solid phases, respectively. The 
thermal and mass diffusion coefficients as, a>£, Dg, DJ are assumed constant. The 
following conditions are usually assumed throughout; 
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In many cases, it is also assumed 

(h) y(t) = a t . 

Equation (a) describes the condition that the original mixture is all liquid at 
temperature T© and concentration Cq. Equation (b) is a consequence of the semi- 
infinite nature of the mixture so that at any time t, the portion near infinity is un- 
changed. Equation (c) assumes that at the solid-liquid interface plane there is an 
interface temperature Ti(t) and that both the solid and liquid phases at x = y(t) have 
this temperature. There is no discontinuity in temperature. Equations (d) and (e) 
state that the concentrations of solid and liquid at the interface are given by the solidus 
and liquidus curves, respectively, of the constitutional diagram for the alloy. Equa- 
tion (f) connects the derivative of the moving boundary with the redistribution of temper- 
ature and Eq. (g) connects the same boundary with that of concentration. Equation (h) 
relates the position of the interface boundary to the solidification time t. 

The conditions on Tg(0,t) and Cg(0, t) are not fixed in our discussion, and a 
number of alternatives are considered; 

1. Tg(0,t) = T^(t) with Tj^(t) equal to a constant for all t; 

2. linear, T^(t) = + t(T^ - T^) /s for t < s and T^(t) = T^ for t > s; 

/ Q 

3. exponential, T^(t) = T^^ + (T^ - T^) e ' so T^(0) = T^ and T^ (w) = T^. 

For Cg(0,t) the conditions considered are C (0,t) = C- usually taken C (T„) or 

at times a condition conserving mass between 0 and oo. s 2 

The two approaches we have pursued may be designated as analytic and numeri- 
cal. The numerical approach can be applied to all three conditions on temperature 
whereas the analytic approach holds oidy the case of constant temperature instan- 
taneously applied, A variant of this analytic method to apply to linear varying tem- 
perature has been investigated. 

An analytic solution to the coupled partial differential equations (1) and (2) sub- 
ject to the initial and boundary conditions (a) through (g) has been given (Ref. B6). A 
numerical program has been designed for the analytic solution. 

These numerical programs developed under NAS 8-27891 are based upon finite 
difference approximations of the partial and ordinary derivatives and involve a variable 
spacing (for improved computing efficiency). The programs have given acceptable 
results and compared well with the reference analytic solution, where comparable. 

The basic physical properties such as densities, diffusivities, specific heats, thermal 
conductivities, and heat of fusion have been held to be constant, and independent of 
temperatures and concentrations. 
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APPENDIX C 


IMPROVED COMPUTER PROGRAM 


The complete computer program for the generalized solidification problem is 
listed herein, together with a glossary explaining the special names used in the 
program. This computer program has the following imique features: 

• Surface evaporation, and its related effects such as material loss, evapora- 
tive segregation, and surface cooling due to the heat of evaporation, have 
been considered 

• Material parameters such as solid and liquid densities, specific heats, 
thermal conductivities, mass diffusivities, and latent heat of fusion or 
evaporation, are allowed to vary with the temperature and composition 

• Realistic phase diagrams involving curved liquids and solidus lines are 
used 

• Two moving boundaries are involved, i. e. , the evaporative boundary and 
freezing boundary 

• Surface temperature is determined by the combined effect of heat radiation, 
evaporative cooling, and thermal diffusion 

Use of Computer Program 

The computer program works well if the following three input program param- 
eters are properly chosen: 1) time step size (DELT), 2) grid spacing (DELX), and 
3) maximum iteration count (NIT). 

The solidification boundary is sensitive to the grid spacing. This is because in 
passing through a mesh point, discontinuity in the computation occurs for the following 
reasons. We compute the derivatives in terms of the temperatures and concentrations 
at the discrete mesh points. When one mesh point is dropped because solidification 
occurs near it, the derivative based on a substituting new mesh point is discontinuous 
with that based on the previous mesh point. Though this discontinuity can be reduced 
by using a smaller time step size, it would be a self-defeating strategy. An alterna- 
tive is to accept the discontinuous results as they occur, advantages being taken of the 
fact that the program corrects itself. Mthough the derivative dz/dt is large when the 
solid-liquid interface passes through a mesh point, it becomes smaller thereafter 
thereby correcting the solidification boundary position. 

The frequency of this self-correction depends on the grid spacing. Too small a 
grid spacing would cause too frequent self-corrections. Too large a grid spacing, on 
the other hand, would obscure the rapid temperature variations around the solidifica- 
tion boundary. This indicates that a proper choice of the grid spacing is required to 
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achieve an optimal tradeoff between accuracy and computing time. There is another 
tradeoff between the time step size and maximum iteration count for optimal comput- 
ing results (see Tables C1-C3). 

Since each evaporation- solidification problem represents a different and unique 
physical situation, each case must be dealt with separately. However, based on our 
experience, the following guidelines would be helpful: 

The first consideration for the choice of the grid spacing is the behavior of the 
evaporation boundary after solidification begins. If the evaporation boundary is vir- 
tually stationary as compared to the solidification boundary, the grid spacing should 
be chosen so that the evaporation boimdary is within the first mesh interval (between 
the first and second mesh points). If, on the other hand, the evaporation boundary is 
moving at velocities comparable to those of the solidification boundary, then the grid 
spacing can be selected more freely. The major consideration in this case is the re- 
lationship behveen the grid spacing and the time step size. For a fixed time step size, 
the grid spacing should be chosen so that at least four time intervals (of step sizes) 
occur before the solidification boundary passes through a mesh point. 

In cases where the evaporation boundary is virtually stationary, one must ex- 
periment to determine an optimum time step size in terms of accuracy and computing 
time. The conditions of the experiment are as follows. Set both the minimum time 
step size (DELTM) and the time printing interval (DELP) to zero. Setting the time 
printing interval to zero will cause the computer to print out every computer time step. 
Setting the minimum time step size to zero will not cause the program to cut back in- 
definitely but will use, as the minimum, the time step size divided by 1024. By ex- 
amining the computed results, one can see at what time step sizes the program is 
running. By examining the actual iteration count (IT), one can see if the program is 
converging or not. If not converging repeatedly, a smaller time step size is indicated. 
If the program is converging most of the time, then the minimum time step size can 
be set at the level of the most frequent time step size and the actual iteration count re- 
examined to see if the program still converges most of the time. For long runs, the 
time printing interval must not be zero or small, but must be chosen in consideration 
with the amount of the required output. 

To impro\'e the computing time on long runs, one should consider enlarging the 
grid spacing as suggested above as one of the tradeoffs. In addition, one may cliange 
the maximum iteration count upwards or downwards to also improve the computing 
time. 


Our computer program has the capability for assuming equal or unequal (doubling) 
mesh point spacings. Our experience, as indicated in Tables C-1 to C-4, shows that 
the unequal spacing scheme gives practically the same accuracy with far less compu- 
tations as compared with the equal spacing scheme. This may be due to the rapidit}^ 
at which the temperature declines at the evaporation boundary. Other physical situ- 
ations may give different results and may indicate that the equal spacing scheme should 
be used. 

The program input parameters consist of a set of integers IX, 1AM, NIT, IM, and 
NONCN; and a set of real numbers DELX, CELT, DELTM, DELP, TF, and S. IX-is 
the maximum number of mesh points to be used in the program. Present, DC 28. 


lAM is the spacing option indicator. If lAM equals 0, the points of mesh are equally 
spaced with grid spacing DELX. If IA.M = 1, an unequal spacing is indicated. The 
first two intervals are equal and set to DELX. Thereafter, each interval is double 
the previous interval in spacing. NIT is the maximum number of iterations as in- 
terpreted in the context of halving the time step size. If the step is begun at the mini- 
mum time step, the NIT is the maximum number of iterations allowed. IM is the 
number of mesh points in actual use. The input value of IM introduces the minimum 
number of mesh points to be used. Thereafter additional mesh points are added as re- 
quired by a substantial change in temperature at next to last mesh point, that is, 1 
degree below the initial temperature, IM is increased until IM is equal to DC. 

Table C-1. Variation of Temperature (*^C) at Evaporative Boundary 
For Various Computation Schemes 


Scheme 

I 

II 

ni 

IV 

V 

Grid 

0.01 

0.01 

0.001 

0.001 

0. 0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 






0. 2 

966.5 

966.5 

966.5 

966. 5 

966.5 

0. 6 

959.4 

959.4 

959.4 

959.4 

959.5 

1.4 

945.7 

945.7 

945. 7 

945.7 

945.8 

1. 8 

938. 9 

938.9 

938. 9 

938.9 

939. 0 

2. 0 

935.5 

935.5 

935.5 

935.5 

935. 6 

2. 05 

934.7 

934.7 

934.7 

934.7 

N. C. 

2.075 

934. 2 

934.2 

934. 2 

934.2 

N. C. 

2. 0875 

934. 0 

934.0 

934. 0 

934. 0 

N. C. 

2.09375 

933. 9 

933. 9 

933. 9 

933. 9 

N.G. 

2.1 

933. 8 

933.8 

933.8 

933. 8 

933.9 

2. 1125 

933. 6 

933.6 

933. 6 

933.6 

933.7* 

2. 1375 

933. 2 

933. 2 

933. 2 

933. 2 

933.3* 

2. 1875 

932. 4 

932.4 

932.4 

932. 4* 

932. 5* 

2. 2875 

930. 8 

930.8 

930. 8 

930. 8* 

930. 9* 

2.4875 

927, 7 

927.7 

927.2 

927.7* 

927.7* 

2.8875 

921. 5 

921.5 

921.5 

921. 5* 

921, 5* 

3. 6875 

909.3 

909.3 

909.3 

909.3* 

909.3* 

5. 2875 

885.8 

885.8 

885.8 




*Hand interpolations 
N. C. not computed 
















Table C-2. Variation of Position ( m) of Evaporative Boundary 
With Time For Various Computation Schemes 


Scheme 

I 

II 

III 

IV 

V 

Grid 

0.01 

0.01 

0.001 

0. 001 

0. 0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 






0.2 

0.122 

0.122 

0.122 

0.122 

0. 12.2 

0. 6 

0.351 

0.351 

0.351 

0. 351 

0.351 

1.4 

0.752 

0.752 

0.752 

0.752 

0.752 

1.8 

0.927 

0.927 

0.927 

0. 927 

0. 928 

2.0 

1.009 

1. 009 

1.009 

1. 009 

1. 010 

2.05 

1.029 

1.029 

1.029 

1. 029 

N. C. 

2.075 

1.039 

1.039 

1.039 

1.039 

N. C. 

2. 0875 

1.044 

1.044 

1.044 

1. 044 

N. C. 

2.09375 

1.046 

1.046 

1.047 

1.047 

N. C. 

2. 1 

1.047 

1.047 

1.047 

1.047 

1.050 

2. 1875 

1.048 

1.048 

1.048 

1. 048* 

1. 052 

2.2875 

1.049 

1.049 

1.049 

1. 048* 

1.053 

2. 4875 

1.051 

1.051 

1.051 

1. 050* 

1. 055 

2.8875 

1. 055 

1.055 

1.055 

1. 054* 

1.059 

3. 6875 

1.061 

1.061 

1.062 

1.061 

1. 066 

5. 2875 

1.072 

1. 072 

1.072 

— 

— 

^'Hand interpolations 
N. C. not computed 
















Table C-3. Variation of Position (/um) of Solid-Liquid Interface 


Scheme 

I 

II 

III 

IV 

V 

Grid 

Spacing 

0.01 

equal 

0.01 

imequal 

0.001 

equal 

0.001 

unequal 

0. 0001 cm 
unequal 

time, ms 






0. 21 



0.109 

0. 109 

0. 106 

0. 24875 

mgm 


0.179 

0. 180 

0. 408 

0. 28875 

rnKm 

0,229 

0.350 

0.350* 

1. 03 

0.36875 

0. 283 

0.283 

0.866 

0. 860* 

2. 76 

0. 52875 

0. 429 

0.429 

2.291 

N. C. 

N, C. 


Table C-4. Variation cf Temperature (°C) At Solid-Liquid Interface 


Scheme 

I 

H 

in 

IV 

V 

Grid 

Spacing 

0.01 

equal 

0. 01 
unequal 

0.001 

equal 

0.001 

unequal 

0, 0001 cm 
unequal 

time, ms 






0. 21 

933.8 

933.8 

933.8 

933.8 

933. 9 

0. 24875 

927.7 

927.7 

927.7 

927. 7* 

927. 9* 

0. 28875 

924. 5 

921.4 

921.5 

921.4* 

921. 6* 

0. 36875 

909.3 

909.3 

909. 3 

909. 3* 

909. 4* 

0. 52875 

885.8 

885.8 

885.8 



*Ifeiid interpolations 
N. C. not computed 
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NONCN is a nonconvergence option. Failure to converge occasionally is not 
necessarily an indication of unacceptable results. Therefore, it is desirable to con- 
tinue computations and examine the results to see if they are acceptable. This is done 
by setting NONCN to 1. If NONCN is set at 0, the nonconvergent results are not 
printed unless called for by the print interval. If NONCN is -1, the program stops on 
nonconvergent results. 

The quantity DELX is the grid spacing. Equal spacing and unequal double 
spacing both make use of this quantity as indicated in the discussion of lAM. The 
quantity BELT is the maximum time interval (step size) for computation. The quan- 
tity DELTM is the input minimum time interval. The program uses as its actual 
minimum the larger of the quantities DELT/1024 and DELTM. Thus, even if DELTM 
is set at 0, the number of halving on cutting back the time interval is at most 10 
^2l0 = 1024). The program in its presolidification phase starts with its actual time 
step BELTS set to DELT/8 and allows it to build up to BELT by quick convergence. 

On the other hand, near the beginning of solidification, BELTS is allowed to cut 
back to as small as DELT/256 in order to find an acceptable start of solidification. 
After solidification has begun, then the restriction of BELTS is between BELTK and 
BELT. If halving reduces BELTS below BELTK, it is set to BELTK. The quantity 
BELP is the print interval. If BELP = 0, then every time step is printed. TPE is 
the time for outputing results. TPR is set originally to BELP and after printout is 
reset to TPR + BELP. The program prints results if the time TIMEl at the end of the 
time step equals or exceeds TPR. The program does not attempt to set BELTS so 
that TPR = TIMEl. This is only a slight inconvenience when the print interval is 
large as compared to DELTM. Generally, we would like DELTM to be set close to the 
most frequent BELTS provided that failure to converge does not ensue on a regular 
basis. TF is the final time of program. This nieans that if TIMEl equals or exceeds 
TF, no additional time steps are taken. 

The decimal quantity S between 0 and 0. 5 is used to determine closeness to a 
mesh point. If the boimdary point (either evaporation or solidification) is such that it 
exceeds the point that divided the mesh interval surrounding the boundary point in the 
ratio (1-S)/S, then the mesh reference point for computation is moved to the next mesh 
point. The introduction of S is to make the transition due to passing a mesh point less 
abruptly discontinuous. The best values of S are between 0. 05 and 0. 15. For com- 
putations on the solid side of the solidification boundary, we continue to use the old 
mesh points until the boundary point passes the point that divides the new mesh interval 
about the solidification point in the ratio S/(l-S). This strategy causes a gradual 
transition from one mesh point to another. The integers III, II2 are used as reference 
point indicators for the solid and liquid sides, respectively. For the evaporation 
boimdary, 113 is used to indicate which points are used. II4 is used only to indicate 
the first mesh point to the left of the evaporation boundary. 

Typical Computer Input 

The definitions of the various inputs fed into the computer are given in the 
Glossary of Program Parameters. Typical input values are as follows: 



IX 


= 28 = maximum number of mesh points 
lAM = 1, unequal, doubled grid spacing 

NIT = 20, maximum number of iterations 

IM =16, actual number of points in mesh 


NONCN 

= 0, allowing the program to continue when nonconvergence 
occurs with no special printout of these results. 

The alloy phase diagrams are determined from the five constants ET, EA, EB, 
EC, and ED, which define the liquidus Ch and solidus lines C as two functions of the 
temperature, T; 

Cj, (T) 

= ED X (ET - T)^ + EC X (ET - T) 

C,(T) 

= EB X (ET - T)^ + EA X (ET - T) 

In our example of 10 mole percent (Cq = 0. 10) of antimony in germanium initially 
imiform at 970®C (Tq = 970) 

ET 

= melting point of pure germanium = 956° C 

EA 

= 0. 128812 X 10~^ 

EB 

= -0. 82218 X 10“*^ 

EC 

= 0. 466678 X lO"^ 

ED 

= -0. 60466 X 10“^ 

The evaporation constants for the solvent and solute as defined previously under 
"The Equation at the Evaporative Boundary" are; 

AU 

= A =0. 1115 X 10^ 
u 

BU 

= By = 0. 863 X 10^ 

EMU 

= My = 0. 2435 X 10^ 

AV 

= Ay = 0. 1171 X 10^ 

BV 

= By = 0. 1803 X 10^ 

EMV 

= My = 0. 726 X 10^ 

EK 

= K = 5. 833 X 10"^ 
e 
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The diffusion coefficient of the solute in solvent in the solid and liquid states are, 
respectively 

DS = D = 0. 10 X 10“® 
s 

DL = 1^ = 0. 10 X 10"^ 

The density p , and latent heat of evaporation, 7 , of the pure solvent are, re- 
spectively, 

RHO = p^= 5.32 
GAMMA = 7 = 160 

V 

Corresponding values for pure solute are: 

RHOU = "= 6* 68 

GAMMAU = 7 = 39 

u ’• 

The above give two derived quantities: 

ALS = aj = /p^c 

ASS =s? = k /fi c 

S s V 

where 

CEE = c = 0. 740 X 10 ^ = specific heat 

The two input parameters in the surface radiation terms are: 

EE = € = 0. 55 = emissivity coefficient, and 

SIG = CT = 0. 136x 10“7=: stefan-Boltzmaim constant. 

Computer Output 

The first line of computer outputs gives the program input parameters IX, LAM, 
NIT, IM, and NONCN, which are defined previously and also in the "Glossary. ” The 
next two lines of computer output give the phase diagram constants (ET, EA, EB,EC, 
and ED) and evaporation constants (AU, BU, EMU, AV, BV, EMV, and EK), respec- 
tively. The next printouts are for CEE, DS, DL, TO, CO, XKL, RHO, GAMMA, 
RHOU, GAMMAU, EE, SIG, T2, and COO, where T2 is the temperature when solidifi- 
cation begins for the melt of initial solute concentration CO. 
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The computed numbers are then outputed as follows; 

IT = actual iteration count 

IM = number of points in mesh 

111 = grid point reference for solid side of mesh 

112 = grid point reference for ll<juid side of mesh 

113 = grid point reference for evaporation boundary 

114 = grid point reference for point lifter evaporation boundary 

These printouts are then followed by the computed values associated with the 
evaporation boundary: time, location y, concentration C, temperature T, extent of 
points X(IM), current time interval BELTS, dy/dt DYDTl, dC/dt DCDTl, dT/dt 
DTDTl. If solidification has not begun, then there is no information about the solidi- 
fication boundary. Otherwise, we have position of the solidification boundary z com- 
puter language (ZIl), solid solute concentration C (CSl), liquid solute concentration 
Cjt (CLl), temperature T (Til), and rate of movement dz/dt (DZDTl). All decimal 
outputs are five per line with excess going to the next lines. 

Representative Computed Results 

The study of the effect of varying the grid spacing DELX on the computed results- 
is summarized in Tables 1 through 4. The five cases considered are: 

Case I: DELX = 0. 01 cm with equal spacing 

Case II; DELX = 0. 01 cm with unequal spacing 

Case III: DELX = 0. 001 cm with equal spacing 

Case IV; DELX = 0. 001 cm with unequal spacing 

Case V; DELX = 0. 0001 cm with unequal spacing 
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APPENDIX D 


AN IMPROVED COMPUTER PROGRAM 



FILE: BINCH6 FCBTRAN P1 


GRUMMAN 


DATA 


S y S I E ^ 


DIMENSION X(28) ,T(20) ,TT(28) ,C(28) ,CC(28) , 

*TEH(10) ,A(84) 

NAMELIST /IN VAR/CEE, DS, DL ,T0 , CO , X KL , R HO, G AMMA , RHOU , G A MM A U , EE, SIG 
♦,ET,EA,ED,EC,ED,AU,BU,EMU,AV, EV,EMV,EK,IX 
D2 (X,F,Y,G,Z,H) = ((H-G) / (Z-y) - (F-G> / (X-Y) ) / (Z-X) *2. 

ABS1 (X) -=AMAX1 (1 . , ABS(X) ) 

UE (V) =EK* (10.** (AU-BU/( 27 3. 12+V) n /SQRT (EHU*( 273. 12 + V) ) 

VE (V) =E’K*(1^.** (AV-BV/(273.12 + V) ) )/SQRT (EMV*{273. 12 + V) ) 
f S (V) = (EB* (ET-V) +EA) * (ET-V) 

FL (V)= (ED* (ET-V)+£C)*(ET-V) 

XCL(V) =ET-2. *V/(EC + SQRT ( (EC) *»2+U .* (ED)* (V) ) ) 

XCSL (V)=ET*2.* V/ ( ( EC-EA) <-SORT { (EA-EG) **2*4.» (EB-ED) *V) ) 

DFL(V) =- (2.*ED* (ET-V) +EC) 

11 = 1 
10=8 

READ(II,100) IX,IAM,NIT,IM,NONCN 

100 FORMAT (1615) 

N ITH=NIT/2 
NITQ=NITH/2 

N ITL=NITH+NITQ 

64 READdI, 101) ET,EA,EB,EC ,ED 

AQUAN = - (EA-EC) **2/ (4.* ( EB-ED) ) 

101 FORMAT (7E''0,0) 

READdI, ri) AI],DU, EMU, AV,EV,EMV,EK 

READ (II, I'^l) CEE,DS,DL,TO ,C0, XKL,RHO,GAMMA , RHOU , G AMM AU ,E£ , SIG 

ALS=XKL/ (RHO*CEE) 

XKS=1. 1*XKL 
ASS = 1 . 1*ALS/1 . 03 
AS=SQPT (ASS) 

AL=SQRT (ALS) 

- READ(II,101) DELX,DELT,DEITM,DE1F,TF,S 
DE1TK=AMAX1 (DELTM,DELT/1024.) 

T2=XCL (CC ) 

COO^FL (T2) 

20 1 DO 1 1 = 1 ,IX 
IF(I-2) 2,3,4 

2 X(1)=0. 

GO TO 1 

3 X (2) =DELX 
GO TO 1 

4 IF(IAM) 5,5,6 
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5 X (I) =X (1-1) +DELX 
GO TO 1 

6 X (I) =X (1-1) +X (1-1) 

1 CONTINUE 

999 HRITE (10, 100 ) IX ,IA M , NIT ,I M, NONCN 
WRITE (10, 102) ET,EA ,EB,EC, ED 
WRITE (10, 102) AU,BU, EMU, AV,BV,IMV 
WRITE (10, 102) CEE,DS,DL, T0,C0,XK1 
*C0^ 

TPR=DELP 

RAT=1. 

TS11 = T0 

TI1=T0 

CSL1=C0 
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FILE; BIHCH6 


FORTRAN PI 


GRUMMAN 


DATA 


SYSTEM 


Y1 1='' . 

ZI1=''. 

II1 =2 
112=2 
113 = 2 
D2T2=0. 

D2T3=T. 

D2T4=0 . 

D2C3=0 . 

D2Ctt=0. 

D2C2='^. 

D2T1=T. 

D2C1='^. 

DTLDX=0. 

DTSDX=0. 

D2C5 = '^. 

D2T5 = '>. 

D2C6 = 0. 

D2T6=n. 

TIME =0. 

DCDX0='^, 

DTDX0=0. 

DELTS=DELT/8. 

TIME1=TIME-»-DELTS 
DO 10 1= 1,IM 
C(I)=C0 
CC (I)=C0 
TT(I) =T0 

10 T(I)=TO 
IFL=0 
IFS=0 
111=112 

Itt IT=P 

IF (IFS.EQ, 1) GO TO 1 99 
IF(TFL) 11,11,20 

11 UO = UF(TO) 
vn=VE (TP) 

IF 1-1 

199 IF(IFn.NE.n) CSL1=CS1 
IF (IFS.EO. 1) IFS = 2 

DYDT<1=U0*EMU*CSLl/RHQU't- V"-*EMV* (1. -CSL1) /RHO 
HB0=-ES*‘SIG*( 273. 1 2+TSII ) ♦ = 4-00’*'GAMMAU*CSL1-VO*GAMMA* (1.-CSL1) 
DCDTO = DCDXO*DYDTP- (UD-VO) *CSL 1 
DTDT0 = DTDX''’*'DYDT j + HBO 
20 YI2 = YI1 + DELTS*DYDT'; 

IF (IFS.EQ.O) ZI2=YI2 

IF (IFS .NE. 0) ZI2=ZI1 +DELTS* DZDTO 

IF(ZI2.GT. X (II 2+1) ) ZI2= (X (II 2) +YI2) /2. 

TSI2 = TSI1 + DELTS+DTDT') 

TI2=T5I2 

CSL2=CSL1 +DELTS*DCDT- 
IF(IFS.EQ.O) CL2=CSL2 
77 IF (IFS. EO.O) GO TO 777 
IF (IFL.GT.I)GO TO 877 
IIT=9 
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DATA 
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877 CS2-FS(TI2) 

CL2 = PL (T12) 

777 D2C2=D2 (ZI2,CL2,X (II2) , CC (II2) ,X(II2*1) ,CC(II2+1)) 

CC(II2) = (CC (II2)+C(II2) + .5*DELTS*(D2C1 + D2C2) *DL ) /2 . 
CALL MOTON(X (112) ,CC(II2) , ZI2 ,C12 , X (112'^ 1 ) ,CC(II2+1)) 

IP (ZI 2 .lt. X(IT2-1) ) CC(II2-1) =CC(II2) ♦(X(II2-1)-X(II2n ♦ 

♦ (CL2-CC(II2) )/(ZI2-X(I12) ) 

IF (II1-II3-1) 83,87, 84 

8T XP=YI2 
CP=CSL2 
GO TO 184 
84 XP = X(II1-2) 

CP=CC(II1-2) 

184 D2C4=D2 (XP ,CP, X(irl-1) ,CC(II1-1) ,ZI2,CS2) 

CC (II1-1) =C(II1-1) + .5*DELIS*(D2C3*C2C4)*DS 
CALL MOTON (X (1 11-1) ,CC(II1-1) , XP, CP ,ZI2, CS2) 

83 IF(II2.EQ.II1.0R.ZI2.LT.K(II2-1)) GC TO 85 
IF (II2-II3.GT. 1) GO TO 185 
XP=YI2 
CP=CSL2 
GO TO 51 

185 XP=X(Il2-2) 

CP=CC(II2-2) 

51 CC(II2-1) =CP + (X(II2-1) -XP)*(CS2 -CP)/ 

♦ (ZT2 -XP) 

8 5 D2T2=D2(ZI2,TI2,X(II2) ,TT (II 2) ,X(II2+1) ,TT (112+1) ) 

IF (D2T2.GT. 0.) D2T2=0. 

TT (II2)= (TT (II2) +T (II2) +.5*DELTS* (E2T 1 + D2T2) *ALS ) /2. 
CALL MOTON (X(II2) ,TT (II2) ,ZI2 , TI2 , X (II2 + 1 ) ,TT(II2 + 1)) 

IF (IFS.EQ.7)GO TO 485 

485 IF (ZI2.LT. X(II2-1) ) TT(II2-1)=TT (II2) +(X (II2-1 ) -X (II2) ) * 

♦ (TI2-TT (II2) )/ (ZI2-X (II2) ) 

IF(II1-II3-1) 69,169,269 

169 TP=TSI2 
XP= YT2 
GO TO 16 
269 TP=TT(II1-2) 

XP=X(II1-2) 

16 D2T4=D2(XP,TP, X(II1-1) ,TT (II1-1) ,ZI2,TI2) 

TT(II1-1)= T (II1-1) +.5*DELTS* (E213 + C2T4) ♦ASS 
CALL MOTON (X (I1 1-1) ,TT (Ill-I) , XP , TP ,ZI2 , TI2) 

69 1;F(II2.EQ.II1.0R.ZI2.LT.X(II2-1) ) GC TO 86 

52 IF (II2.LT. II3-1) GO TO 186 
XP=YI2 

TP=TSI2 
GO TO 352 

186 XP=X(Il2-2) 

TP=TT (II2-2) 

352 TT(II2-1)=TP + (X (II 2- 1 ) -X P) * ( TI 2 -TP)/ 

*(ZI2 -XP) 

86 IF (IFS.EQ.O) GO TO 299 

DCLDX= (CL2-CC(II2) ) / (ZI2 - X ( II 2) ) 

IF (D2C2.GT.0.) DCLDX=DCLEX-D2C2* (X (II2) -ZI2) /2 . 

IF (D2C2.LT.0.) D2C2=0. 

DTLDX= (TI2-TT (II2) ) / (Z I 2-X ( II 2) ) 
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GRUMMAN 


DATA 
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IF (D2T2.lt. 0.) DTLDX=DTLDX-. 5*D2I2^ (X(II2) -5112) 

IF (II3.EQ.II1) GO TO 386 
XP=X(II1-1) 

TP=TT(II1-1) 

CP=CC(II1-1) 

GO TO 486 
386 XP=YI2 
TP=TSI2 
CP=CSL2 

DCSDX= (CS2-CP) / (ZI2-XP) 

D?SDX= (TI2-TP) /(ZI2-XP) 

GO TO 686 

486 DCSDX= (CP-CS2)/(X'P-ZI2) -D2C4* (XP-ZI2) /2. 

DTSDX= (TP-TI2)/(XP-ZI2) -D2T4* (XP-ZI2) /2. 

686 DZDT= (DL*DCLDX-DS*>DCSDX)/ (CS2-CL2) 

DZDTT=(XKS*DT3DX-XKL*DTLCX)/ (RHC*GAMMA) 

FSL=DZDT* (CS2-CL2) /DZDTT 

IF (FSL.GT.C. .08.FSL.lt. AUUAN ) GO TO 772 

TII=XCSL (FSL) 

GO TO 771 

772 TII=FT 

771 C0E1 = XKS/ (ZI2-XP) 

C0E2=XKL/( X(II2) - ZI2) 

586 TI= (RH0*GAMMA*DZDT + C0E1* (TP-D2T4*. 5* (XP-ZI2) **2) +C0E2* (TT (II2) - 
♦ ♦D2T2* (X (II2)-ZI2)**2) ) /(C0E1+C0E2) 

773 IF (TI.LT.TSI2 .AND.TI.LT.TII) TI=TII 

IF (TI.LT.TSI2 .0R.TI.GT.TT(II2 + 1) )TI=TT(II2) 

IF (DZDT.LT.O. .AND.DZDTT. GT.O.) D2DT=DZDTT 
IF (A8S (TI-TI2) -1. E-5»ABS1 (TI+TI2) ) 58 7,587,770 

587 IF (ABS (DZDT-DZUTI) -1.E-3*ABS1 (DZDT + DZDT1) ) 298,298,770 
770 TI2=(TI + TI2) /2. 

IF (DZDT .LT.O.) DZDT=DZDT 1 
DZDT1= (DZDT+DZDTI) /2. 

ZI2=ZI1+. 5*DELTS* ( DZDT1 +DZDT0) 

IF (ZI2. GT. X (II2-M) ) ZI2= (X (II2) + YI2) /2. 

IIT-IIT+1 

IF (IIT.GT. 5) GO TO 298 
GO TO 877 ' 

298 IF(TI.GT.TT{II2) )TI=TT (II2) 

IFI.TI.lt. TSI2) TI=TSI2 
IF{II1-II3-1) 398, 498,598 

59E D2T6=D2 (II2,TSI2,X (II3) ,TT(II3) ,X (113 + 1) ,TT (II3+1) ) 

D2C6=02 (YI2,CSL2,X (II3) ,CC(II3) ,X(II3+1) ,CC(II3+1)) 

DtDXl= (TSI2-TT (II3) ) /{YI2-X (II3) ) -D2T6* (X (II3) -YI^) /2. 

DCDX1= (CSL2-CC (II3) ) / ( Y I2-X ( 1 13) ) -D2C6* ( X ( II3) -YI2) /2. 

GO TO 599 

498 DTDX1=DTSDX+D2T4* (XP-ZI2) 

DCDX1 = DCSDX+D2C4* (XP-ZI2) 

D2T€=D2T4 
D2C6=D2C4 
GO TO 599 
398 DCDX1=DCSDX 
DTDX1=DTSDX 
GO TO 599 

299 DTDX1= (TT (II2) -TI2) / (X ( II2) -YI2) -. 5*D2T2* (X (II2) - Yl2) 
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DCDX1=(CC(II2) -CL2) /(X(II2) -YI2) 5*D2C2* (X (II2) -Yl2) 

599 U1=UB(TSI2) 

V1=VF(TSI2) 

HB1=-EE*SIG* (273.12+TSI2) ♦♦4-Ul*GAMMAU*CSL2-Vl*GAKMA* (1 . -CSL2) 

DYDT1=U1 *EHU*CSL2/RH0U*V1*EMV* ( 1 . -CSL2 ) /RHO 

DCDT1=DCDX1*DYDT1- ( U1-V1) *0312 

DTDTl=DTDXr*DYDTl+HBl 

IF {IFL. LT. 2 ) GO TO 76 

IF(IFS.GE, 1) GO TO 174 

IF(TSI2.GT.TT(II2) ) GO TO 399 

TT2=XCL (CSL2) 

IF(TSI2.GT.TT2) GO TO 174 
IF (TSI2.GE.TT2-. 05) GO TO 400 
IF (DELTS. LE.DSLTK) GO TO 400 
GO TO 399 

76 YI = YI1+.5* (DYDT0+DYDT1) ♦BELTS 
IF(IFS.EQ.O) ZI = YI 
TSI =TSI1+.5*(DTDTC+DTDT1) *DELIS 
IF (TSI.LT.O.) TSI=TSI1 
CSL=CSL1 + . 5* (DCDTO+dCDTI) ♦DELTS 
IF (IFS.EQ.T) GO TO 73 
IF (IIT.GT.5) GO TO 70 

IF(TSI.GT.TT (II2) .OF.TSI. GT.TI2) GO TO 70 

73 IF (ABS (YI-YI2) -1.E-6*ABS1 (YI + YI2) ) 74,74,70 

74 IF(ABS(CSL-CSL2) -1.E-4*ABS1 (CSL+CSL2) ) 75,75,70 

75 IF (ABS (TSI-TSI2) -1. E-5>^ABS1 (TSI+TSI2) )7,7,70 
70 TSI2=(TSI + TSI2) /2. 

IF (IFS.EQ.O) TI2=TSI2 
YI2= (YI+YI2) /2 . 

IF (IFS.EQ.O) ZI2=YI2 
CSI2= (CSL+CSL2) /2. 

IF (IFS.EQ.O) CL2 = CSL2 

IF(ZI2.LT. (1 .-S) ♦X(II2) +S*X(II2-1) ) GO TO 24 
22 IF(II1-II2) 24,96,96 
9 6 112=112+1 

46 D2C1=D2 (ZI1,CL1 ,X (II2) , C ( II2) , X (II2 + 1 ) ,C (II2+1) ) 

D2T1=D2 (ZI1,TI1,X (II 2) ,T(II2) , X ( II 2+ 1) , T (II2+ 1) ) 

24 IF (IT-NITQ) 48,174,160 
16'' IF(IT-NITH) 48,47,48 

47 IF (DELTS-DELTK) 48,48, 53 
53 DELTS=DELTS/HAT 

152 IF(FAT-1.) 153,153,154 

153 RA1=2.*RAT 

DELTS=DELTS/2. 

GO TO 152 

154 TI«E1=TIME+DELTS 
212 DO 45 I=II3,IH 

TT(I)= (TT(I) +(RAT-1.)*T (I) )/RAT 
4 5 CC (I) = (CC (I) + (RAT-1.)*C(I))/RAT 
148 RAT=1 . 

48 IT=IT+1 

IF (IT.EQ. NITH+ 1) GO TO 20 
IF (IT-NITL) 161,174,161 
161 IF (IT-NIT) 77,77,402 
402 IFCDatS.GT. DELTK) GO TO 399 
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BIN02370 
BIN02380 
B1N02390 
BIN02400 
BIN 024 10 
BIN02420 
BIN02430 
BIN02440 
BIN02450 
B1N02460 
BIN02470 
BIM02460 
BIN02490 
BIN02500 
BIN02510 
BIN02520 
BIN02530 
BIN02540 
BIN025^ 
BIN02560 
BIN02570 
BIM02580 
BIN02590 
BIN02600 
BIN02610 
BIN02620 
BIN02630 
BIN02640 
BIN02b50 
BIN02660 
BIN 02670 
BIN02680 
B1N02690 
BIN02700 
BIN02710 
BXN02720 
BXN02730 
BXN02740 
BIN 02750 
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FILE: BINCB6 FOHTRAN P1 


GRUMMAN DATA SYSTEM 


IPl=2 
GO TO 77 

26 IF (NONCN) 99,U B,48 
7 IFl=IFL-M 
DZDT1=DZDT 

IP (IFS.NE.O) ZI2 = ZI1 +. 5*DE1TS* (DZDTO + DZDT1 ) 

TSI2=TSI 

YI2=YI 

IF (IFS.EQ. 0) ZI2 = YI2 
CSL2=CSL 

IF (IFS.NE.O) TI2=TI 
IF (TI2.LT.TSI2) TI2=TSI2 
IF (IPS. EQ.O) TI2=TSI2 
IF(IFS.EQ.O) CL2=CSL2 
GO TO 77 

399 IF(DELTS.LE.DELTK) GO TO 26 
IT= NITH 
IFL = 1 
GO TO 153 

117 IF (TIME. EQ. TIMED TIME1 = TIME1 ♦DELIS 

TIME=TIHE1 

IF (RIIT, NE.2. . AND. RAT. NE.O.) DELTS=DELTS/E AT 
RAI=1 . 

I?(IT-HITQ) 02,82,81 
82 IF(DBLTS.GT.DELT/2.) GO TO 81 
DSLTS='DELTS + DELTS 
RAT=2. 

81 TIME1=TIHE1+DELTS 
282 D2C1=D2C2 
D2C3=D2C4 
D2C5=D2C6 
D2T3=D2TU 
D2T1=D2T2 
D2T5=D2T6 
DYDT0 = DYDT1 
DCDTO=DCDT 1 
U0 = 01 
V0 = Y1 

DTSDX1=DTSDX 
DCSDX1=DCSDX 
IF(IFS.NE.O) DCDXC---DCDX1 
IF (IFS. NE.O) DTDX'^ = DTDX 1 
IF (IFS.NE .0) DZDT0 = DZDT1 ' 

CSL1=CSL2 

TSI1=TSI2 

YI1=YI2 

ZI1=ZI2 

IF (IFS.NE.O) CS1 = CS2 

HB0=HB1 

C11=CL2 

TI1=TI 2 

IF (YI2.LT. X(II3) ) GO TO 410 
II3=II3+1 

D2T5=D2(YI1,TSI1 ,X(II3) ,TT(II3) ,X(II3+1) ,TT(II3+1) ) 
D2C5=D2(YI1,CSL1,X (II3) ,CC (1I3) ,X(II3+1) ,CC (II3 + 1)) 


BIN02760 

UIN02770 

BIN02780 

BIN02790 

BIN C2800 

BIN02810 

BIN02620 

BIN02830 

BIN02B40 

BIN02850 

B1K02B60 

BIN02870 

BIN02880 

BINU2890 

BIN02900 

BIN02910 

BIN02920 

BIN 02930 

BIN02940 

BIN02950 

BIN 02960 

BINU2970 

BIN02960 

BIN02990 

BIN03000 

BIN03010 

BIN03020 

BIN 030 30 

B1N03040 

BIN0305U 

BIN03060 

BIN03070 

5IN030B0 

BIN 03090 

BIM031C0 

BIN03110 

BIN03120 

BIN03130 

BIN 03 140 

BIN03150 

BIN03160 

BIN03170 

BIN03180 

BIN03190 

BIN03200 

BIN 032 10 

BIN03220 

BIN03230 

BIN 03240 

BINU3250 

BIN03260 

BIN03270 

BIN03280 

BIN03290 

BIN03300 




’’HE; BINCR6 FORTRAN P1 GRUhWAN DATA 

SYSTEM 


• 410 I? (Ill . EQ.II2.OR.ZI2 .lt. S*X (II2) ♦( 1.-S) *X(II2-1) ) GO 10 33 

BIN03310 

■ ^ 

111=112 

BIN 03320 


IP(IH-II3-1) 33,310, 110 

BIN03330 


310 xp=yri 

BIN03340 


CP=CSL1 

BIN 03350 


tp=tsii 

BIN03360 

.. ^ 

GO TO 210 

BINJ3370 


no xp=x(iii- 2 ) 

BIN03380 


CP=CC(II1-2) 

BIN03390 


TP=TT (III -2) 

BIN034C0 


2ir D2C3=n2(XP,CP,X(I11-1> ,CC{II1-1) ,ZI1,CL1) 

BIN03410 


D2t 3=D2 (XP,TP, X(II1-1) ,TT (Ill-I) ,ZI1,TI1) 

UIN03420 

^ * 

IF(II1-II3.NE. 1) GO TO 33 

BIN03430 


D2T5=D2T3 

BIN03440 

. : 

D2C5=D2C3 

BINC3450 


GO TO 33 

BIN034b0 

■'"I 

174 111=112 

U1N03470 


IP (ZI2.lt. X(II2-1) ) III=III-1 

BIN034BO 


IP (II1-II3.LT. 2) GO TO 1 8 

BINC3490 


TP (111 -II 3. GT. 2) GO TO 29 

BIN03500 

- ■? 

TT (II3)= (T(II3) /DELTS+. 5*ASS*D2T5 

BIN03510 

■■ J 

. * + TT(II3+1) ♦RSS/(X(II3+1)-X (II3) )/ (X (II3 + 1) -YI2) +TSI2*ASS/ (X 

HIN03520 


*ai3)-YI2)/(X(II3+1)-YI2) )/(1./DELIS*ASS/(X(II3 + 1)-YI2)*(1./(X(II3BIN03530 


*»1)-X (II3) ) + 1. / (X(II3) -YI2) ) ) 

BIN03540 


CC(T13) =(C (II3) /DELTS+. 5*ES>kD2C5 

BIN03550 


*+CC (II3^1)*DS/ (X (II 3 + 1) -X (II3) ) /(X (II3+1) -YI2) +CSL2*DS/ (X (I 

BIN03560 


*I3)-YI2)/(X(II3 + 1)-YI2) ) / (1 ./DELTS + DS/ (X (II3+1) -YI2) *(1./(X (113 + 1) B1N03570 


*-X(II3)) +1./(X (TI3) -YI2)) ) 

BIN03580 

' T 

GO TO 18 

BIN03590 

5 

29 X(II3-1) =(YI2+YI1) /2. 

UINJ3600 

^ i.- 

C(II3-1)=CSL1 

BIN03610 


CC(II3-1) =CSL2 

BIN03620 

•T " 

T (TT3-1) =TSI1 

BIN03630 

I 

TT(II3-1) =TSI2 

BIN03640 


CALL TRIST(X,T,TT ,113-1 ,111-1 , ASS , BELTS , A) 

aiN03650 


CALL TRIST ( X, C , CC, II3- 1, III- 1 , DS, CELTS, A) 

BIN03660 


GO TO 18 

BIN03670 

V 

27 IF(T(in-1) -T0 + 1.E0) 30,32,32 

BIN03680 


30 IP (IM-IX) 15,32,32 

BIN03690 


TS IM=IM+1 

BIN03700 


T(IH) =TC 

BIN03710 


C(IH)=CO 

BIN03720 

'■ • 

32 TT(IM) =T(IM) 

BIN03730 


CC (IM)=C<IH) 

BIN03740 


GO TO 14 

BIN 03750 


18 IF (BELTS. LT. DBLTK) DELTS = DELTK 

BIN03760 


CALL TRIST(X,T,TT,III,IM,ALS,DELTS, A) 

BIN03770 


CALL TRIST (X,C,CC, III, IH, DL, BELTS, A) 

BIN03780 


DO 21 I=II2,m 

BIN03790 


IF (I.EQ.II2) GO TO 21 

BIN03800 

*• ' 

21 CONTINUE 

BIN03810 


219 IF (IPL.EQ.2) GO TO 117 

fiIN03820 

' i 

GO TO 48 

DIN03830 


33 IF(TIME -TPR) 50,34,34 

BIN03840 

* *' 

5v0 IP (NONCN) 98,98,54 

BIN03850 

: I 
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FILE 


BINCE6 FORTRAN Pi 


G B U tJ M A N 


DATA 


S y S T E « 


5U 

IP(IT-NIT) 98,38,42 

BIN03860 

m 

IF (IIELP.SQ.C.) GO TO 42 

BIN03870 


TPB=TPR+DELP 

BINU5880 

42 

WRITE (10, lO'^) IT, IN ,111,112 ,113 

B1N03890 


WRITE (10, 102) TIME, yi1 ,C SL1 ,TSI1 ,U1 ,Vl ,HB1 ,DYDT1 ,DTDT1 ,0CCT1 ,DELTS 

BIN03900 


IF (IFS.NE.O) WRITE (10, 10 2) Z1 1, CS 1,CL 1,T1 1 , DZDT1 

DIN03910 

143 

WRITE (10, 10 2) (TT(I) ,I = II3,IM) , (CC(I) ,I=IIJ,Xrt) 

BIN03920 

T*? ^ 

F0RMAT(5E14. 6) 

BIN03930 


IFC^IME -TF) 98,99 ,99 

BIN03940 

98 

DO 97 I=II3,IH 

BIN03950 


TTT=TT (I) +RAT* (TT (I) -T ( I) ) 

BIN03960 


T (I)=TT(I) 

BIN03970 


TT(I) =TTT 

BIN03980 

97 

CONTINUE 

BIN03990 


DO 197 I=II3,IM 

BIN04C0O 


CCC=CC (I) +RAT* (CC(I) -C (I) ) 

BIN04010 


C(I) =CC(I) 

BIN04O20 

197 

CC (I) =ccc 

DING 40 30 


IFL = 1 

BIN04040 


GO TO 27 

BINO4050 


CSS=FS (T12) 

BIN04060 


CLL=FL (TI2) 

BIN04070 


DEI7.= (CLL-CI.2)/ ( (CLL+CC (II2) ) -2. *CSS) *■ ( K (II2) - YI2) 

BIN04080 


ZI2=YI2*DELZ 

BIN04090 


DZDTO = DELZ/DELTS 

BINO41C0 


DZDT1=DZDT'^ 

BIN04110 


CL2=CLL 

UIN04120 


CR2 = CSS 

BING4130 


DCDX''=0. 

BIN 34 140 


DTDX0=0. 

BINe4150 


L'TDX1=0. 

BIN04160 


WRITE (10, 103)TIME1, YI2,TI2,ZI2,CS2,CL2, DELZ , DEBTS , DZDIO 

BIN04170 

1''3 

FORHAT(‘ SOLIDIFICATION HAS BEGUN V(5E14.6)) 

BIN04180 


IFS = 1 

BIN04190 


GO TO 174 

BIN04200 

99 

PFADdI, INVAR) 

BIN J4210 


IF (IX.GT.O) GO TO 999 

BIN04220 


STOP 

BING4230 


END 

BIN 04240 



GLOSSARY OF PROGRAM PARAMETERS 


AS, AL 
ASS, AIR 

AU, BU, AV, BV, EK 

CEE 

CEEUO 

CEEUT 

CEE VO 

CEEVT 

C^ 

CLl, CL2 
COEl 
COE2 
COO 

CSl, CS2 
CSL, CSLl, CSL2 

D2C1, D2C2 
D2C3, D2C4 


= temperature diffusion coefficients 
= squares of temperature diffusion coefficients 
= evaporation constants for solute and solvent 
= C specific heat 

= c^^ - constant term for specific heat of solute 
= c^^ - linear coefficient for specific heat of solute 
= c^^ - constant term for specific heat of solvent 
= c^rj, “ linear coefficient for specific heat of solvent 
= concentration of solid at the start of solidification 
= concentration of liquid at solid-liquid boundary 
= k divided by the interface to the first left mesh point 
= divided by the interface to the first right mesh point 

= equals CO 

= concentration of solid at solid-liquid boundary 
= concentration at evaporation boundary 

at liquid side of solid-liquid boundary 
z 

at solid side of solid-liquid boundary 




DIO 


at evaporation boundary 


D2C5, D2C6 
D2T1, D2T2 

D2T3, D2T4 

D2T5, D2T6 
DSC 2 

DSC 4 

DST2 

DST4 

DCDTO, DCDTl 
DCDXO, DCDXl 

DCLDX 

DCLX2 



at liquid side of solid -liquid boundary 



at solid side of solid-liquid boundary 


at evaporation boundary 


= third derivative approx. 


ax® 


= third derivative approx. 


a®T 

ax® 


= third derivative approx. 




= third derivative approx. 


a^T 

ax^ 


dC 


^ derivative of concentration at evaporation 
boundary 


dx/y 


partial derivative of concentration at 
evaporation boundary 


partial derivative of concentration in 
liquid at solid-liquid boundary 

= a c ^ / a X 



Dll 



DCSDX 


partial derivative of concentration in 
solid at solid-liquid boundary 



DCSX2 

= d e j / d X 


DELP 

= time print interval 


DELSC 

= ratio of step size the program wants to use to the 
computed step size 


DELT 

= maximum time interval (step size) 

1 

DELTC 

= indicated step size for the initial amoxmt of solid 

■ s 

DELTK 

= larger of quantities BELT/1024 and BELTM 

■ - 

DELTM 

= minimum time interval 

1 

BELTS 

= current time interval 


BE LX 

= length of the first mesh interval (Xq, x^^) 


BFL (V) 

= derivative of liquidus equation 


BLC 

= B-j-p - linear coefficient for concentration for mass 
diffusivity of liquid 

: 

BLT 

= B, „ linear coefficient for concentration for mass 
diffusivity of liquid 

- ' 

BLO 

= B, ^ - constant term for mass diffusivity of liquid 
LO 

- • 

BS, BL 

= mass (concentration) diffusion coefficient 

- ■, 

BSC 

= Bop - linear coefficient for concentration for mass 
diffusivity of solid 


BSO 

= Bq^ - constant term for mass diffusivity of solid 
bO 

* 

BST 

= Boq, “ linear coefficient for temperature for mass 
diffusivity of solid 

« 

BTBTO, BTBTl 

dT 

= derivative of temperature at evaporation boundary 

: X 
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DTDXO, DTDXl 


DTLDX 

DTSDX 

DYDTO, DYDTl 

DZDT, DZDTT, DZDTO, 
DZDTl 

EE, SIG 

EMU, EMV 

ET, EA, EB, EC, ED 

FS, FL 

G 

GAMMAO 

GAMMAT 

GAMMAU, GMMA 
GAMMUO 

GAMMUT 
HBO, HBl 
lAM 
IFL 

IFS 


partial derivative for temperature at 
evaporation boundary 


partial derivative of temperature in 
liquid at boundary 

partial derivative of temperature in solid 
at boundary 

= ^ derivative of evaporation boundary 

= ^ derivative of solid-liquid boundary 

= radiation constants € , <r 
= molecular wei^t of solute and solvent atoms 
= phase diagram constants 

= arithmetic functions for solidus and liquidus curves 
= temperature gradient, °C/cm 

= Y^q - constant term of lateral heat of fusion of solvent 
= 7 ^ - linear term of lateral heat of fusion of solvent 

= specific heats of solute and solvent 

= 7 - constant term of lateral heat of fusion of solvent 

= 7 uT “ linear term of lateral heat of fusion of solvent 

= heat balance sum of evaporation and radiation terms 

= spacing option: 0 indicates equal, 1 unequal doubling 

= indicator of convergence: 2 on convergence, 

<2 before convergence 

= indicator of beginning of solidification: IFS = 0 before 

solidification 

IFS > 0 after 
solidification 
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ni 

= grid point reference for solid side of mesh 

II2 

= grid point reference for liquid side of mesh 

II3 

= grid point reference for evaporation boundary 

II4 

= grid point reference for point after evaporation 
boundary 

in 

= grid point reference for point after solidification 
boundary 

IM 

= number of points in mesh 

IT 

= actual iteration count 

ITT 

= parameter controlling surface temperature condition 
-1 for instantaneous temperature drop; 0 for linear 
temperature drop; 1 for exponential drop. Ratio 

IX 

= maximum number of points in mesh, < 28 

NIT 

= maximum iteration count 

NITH 

= half of NIT 

NITL 

= 3/4 NIT 

NITQ 

= quarter of NIT 

NONCN 

= nonconvergence option; 1 indicates proceed and 

printout 

0 indicates proceed but 
do not printout 
-1 indicates program stop 

RAT 

= ratio of current time step over previous time step 

RHOO 

= PyQ - constant term of density of solvent 

RHOGAM 

= p y - product of density x latent heat of fusion 
as determined for concentration of solute and 
solvent - liquid 



RHOT 

- li 

RHOUO 

= P. 

uo 

RHOUT 

= ^ rr 

u,T 


f* ^ linear coefficient of density of solvent 


- linear coefficient of density of solvent 


S 


= a measure of the closeness a number between 0 and ^ 
used to determine when computing scheme should 
change the base mesh point. If the ratio of the distance 
from interface (solid -liquid) to the next ri^t mesh- 
point, Xj, to the last mesh interval (xi_]^, x,) is less 
than s, then mesh point xj_i is used instead of xi 
in determining derivatives . 

= temperati-ire of cooler zone 


T_ = temperature at start of solidification 

A 

To = time to freeze ZL = ZL/YD 

O 


T1 

TA 

TF 

TI, TII, Til, TI2 

TIME 

TIMEl 

TIMEO 

TO, CO 

TPR 

TS 

TSI, TSIl, TSI2 
TTF(V) 


= Lower bound for surface temperature variation. 

= ambient temperature 
= final time 

= temperatures at solid -liquid boundary 
= time at beginning of time interval 
= time at end of time interval 
= time of start of solidification 
= initial temperature and concentration distribution 
= time for printing results 

= the time when the lower bound of surface temperature 
is obtained. After this time, the surface temperature 
remains constant at Tl. 

= temperatures at evaporation boundary 

= computed temperature as an exponential function of 
time 
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TU 

= eutectic temperature 


TX 

= proportion of distance between and to that 

between and 

' ! 

UE, VE 

= arithmetic function definition for evaporation rates 

M ^ » 

UO, U1 

= rates of evaporation of solute 

l 

VO, VI 

= rates of evaporation of solvent 


XKLU 

= value of thermal conducti\dty of solute as liquid 

•. ,• 

XKLUO 

~ ]^uo ” term for thermal conductivity of 

solute as liquid 

■ 1 
i' 

XKLUT 

= p - linear coefficient for thermal conductivity of 

solute as liquid 

1 

* i 

XKLV 

= value of thermal conductivity of solvent as liquid 


XKLVO 

~ ^vo ^ term for thermal conductivity of 

solvent as liquid 

- ] 

XKLVT 

= - linear coefficient for thermal conductivity of 

solvent as liquid 

1 

c 

XKS, XKL 

= ks » interphase boundary equation 


XKSU 

= value of thermal conductivity of solute as solid 

- ' 

XKSUO 

= “ constant term for thermal conductivity of 

solute as solid 


XKSV 

= value of thermal conductivity of solvent as solid 

: i 

XKSVO 

= k - constant term for thermal conductivity of 

solvent as solid 

■■ ^ 

YI, YIl, YI2 

= values of y (evaporation boundary) 


YD 

= is the velocity of solidification, i. e. , the velocity of 



the solid-liquid interface. 
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YG 

ZIl, ZI2 
ZL 


= first approximation for the amount of solid 
= values of z (solid-liquid boundary) 

= graded temperature zone length, cm 
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APPENDIX E 


ECONOTECHNICAL MODEL OF CRYSTAL GROWTH 


A. Introduction 


An econo technical model of crystal growth can be set up to study the effects of 
such input variables as growth technique or equipment, alloy system, initial solute con- 
centration, initial melt temperature, and material thermophysical parameters on such 
output variables as solidification velocit 5 q growth morphology, cooling time to solidi- 
fication, evaporative material loss at anytime, crystal size and perfection, power con- 
sumption, and cost. This model may be in the form of a series of linear equations each 
expressing an output variable as a linear combination of the value of the input variables 
multiplied by its unique effect. The model can be more complicated and involves not 
only main effects but also higher-degree terms and multi-factor interactions. The 
values of these main effects, higher-degree terms, and interactions can be determined 
from actual experimental or computer simulation data. Statistical designs are useful 
in planning the experimental or computer runs so that such determinations may be done 
with maximum confidence even with limited available data. 

When fully developed, these correlation equations will allow the prediction of 
crystal size, quality, and cost; the systematic selection of best growth equipment or 
alloy system; the optimization of growth or material parameters; and the maximization 
of zero-gravity effects. 

B. The Model 

Specifically, for small changes A x^ in any of the n (critical or interested) input 

variables xj, any output variable Sj^ can be considered as being linearly related to the 

various values of the x| s orAxjs, i. e. , 

n * _ n 

= 2 a.x. = a X + a.x. +. . .+a x = S, 2 a. Ax. (1) 

. „ iiooii nnk.. 11 ^ ^ 

1=0 1=1 

where the aj's are the respective effects of the n input variables, and Si^ = 3-o^o 
average value of the output variable Sj^ witliin the experimental domain studied. 

An example of use of the models to predict the cooling time to solidification and 
evaporative material loss at beginning of solidification based on initial melt composition 
and temperature will be given later in this section. 
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It can be seen that any effect aj[ for the input variable is the partial derivative 
of with respect to x^, or 



( 2 ) 


If we chose a particular set of conditions (i. e. , values of x.'s) as the standard 
or base line conditions and refer to these x. , *s to determine the A x.'s, then 
Eq. 1 becomes 


Si = S, , 

k k,b^ 


n 

+ 2 

i= 1 




= s, , + 


n 

2 

i = l 




dXi 


Ax, 


(3) 


where S, . is the value of S, at the standard or base level conditions under which all 
. hj k 

the A x.'s are zero. 

1 

In the more general case, particularly when the A x.'s are large, the same out- 
put variable Sj^ may be also related to higher-degree terms of one or more of the in- 
put variables x.. Hence the following improved model for Sj^: 



Further, the input variables in crystal growth are known to interact strongly. 
Thus, there is no best growth equipment or alloy system. But there is a best combina- 
tion of growth equipment and alloy system. Even three or more variable interactions 
are usually present. The classical approach in which only one variable is changed at 
a time, therefore, seldom reaches the true or global optimum conditions and often is 
of limited value. 

Also, the even more general model of crystal growth is; 


II 




= s 


L , - + 2 a.Ax +2b A (X +2 d. A (x.^) + 2 f.. A (x.x.) 
k, b./ 11 1 ' 1 '' 1 ' 1 ^ 1] 1 ]' 


The value of as observed, i, e. , s^, in actual experiments or computer 
simulations, is always different from the true, by measurement or computational 
errors which are assumed to be randomly distributed. Hence, 

Sk,b^= ®k,by + 

and sf = S, , . + 2 a. Ax. + e for the linear model, 
k k,bi» 11 

Let us take the special case when the standard or base level conditions are ground 
level or one-g conditions, while the test conditions are the space or zero-g conditions. 
The zero-gravity effect on the kth output variable is, then: 
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(8) 


ZGE, =As, =Sa.Ax. + (e^^-e)=2aAx + e» 

t k, Og k, Ig ’ ’ ^ ue ' 


where e' = ^2 e = 


%ji 


1 1 
= e. 


hJ 


The procedures for determining these various effects a’s, b’s, d's, fs, 
and g's by least-square techniques from the experimental or computer data are fairly 
straight forward. See, for example. References El and E2. Techniques are also avail- 
able to sequentially determine these same effects or coefficients and automatically 
adding a new important term or eliminating an unimportant term in the model 
equations. 

These model equations are very similar to those used in statistical design of 
modern experiments and the subsequent variance analyses (Refs. El and E2). In fact, 
the actual or computer experimental runs should fully use such modem statistical 
techniques to maximize the yield of reliable information. 

C. Optimization Techniques 


Once these correlation or prediction equations are established, one can estimate 
any output variable S^, with a known degree of confidence, from the known values of 
the input variables. The procedure to optimize a particular Sj^ relative to one or two 
input variables involves simply equating the partial derivative of with respect to 
these variables to zero and checking to see if the extreme S\^ thus obtained is indeed 
the optimum sought. If several or all the input variables are allowed to vary, one 
can even use the method of steepest ascent (Ref. E3) to simultaneously vary the 
input variables x^^'s, to minimize the number of required experimental runs, for 
example, in the expensive space environment. 

We can even simultaneously optimize a number of output valuables S ’s represent- 
ing , e. g. , crystal size, crystal perfection, growth velocity, power consumption, 
contamination index, and cost. These optimizations may be done with or without 
certain constraints on the available ranges of the input variables x^'s, with the aid of 
linear or dynamic programming techniques. 

D. Uses of the Model 

in summary, when fully developed, this program should enable us to; 

1. predict any output variable or growth result upon knowing the values of input 
variables; 

2. detect, analyze, and determine any Interactions among the various input 
variables; 

3. Rigorously and quantitatively compare hvo or more growth methods, equip- 
ments, alloy systems, . . . for a desired output variable; 

4. Define the best or most promising regions of space crystal growth study; 

6. Help to achieve best crystal growth results in space by starting in an optimal 
growth region and zero in to optimization b 5 ^ the method of steepest ascent; 
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6. Systematically select the unique combinations of growth pt^rameters (or 
input variables) for maximized zero-gravity effects; 

7. Finding the most cost-effective approach to crystal growth in space; and 

8. Putting space crystal growth work as a whole on a more objective, efficient, 
and scientific basis. 


E. Example of Usage 


As an example, we used the computer program previously discussed (in Section 
II for calculating the temperature and concentration distributions in a semi-infinite, 
Ge-Sb binary alloy liquid under unidirectional cooling through surface evaporation and 
radiation) to determine the time in seconds to initiate solidification tj^g, and the 
length in centimeters evaporated away before solidification starts, x^^g. We assume 
here that tji^g and X|^g are the output data we want to optimize through controlled varia- 
tions of two other controllable parameters, 1, e. , initial melt temperature in 
and initial melt composition c^ in moles. The use of the techniques described in 


this section, and even of the particular prediction equations for t^^ 
shown. 


and x^g, will be 


The variations of the time to solidification and length of evaporated material 
with T^and c^ are shown below for the system of antimony in germanium: 


_3 

Table E-1 Variation of tj^^ (x 10 sec) with T^ and c^ 


c /T 
o o 

970 

980 

990®C 

Average 

0.01 

0. 9 

1.5 

1.9 

1.4 

0.02 

1. 0 

1.6 

2.1 

1. 6 

0.03 

1.2 

1.7 

2.2 

1. 7 

0. 04 

1.3 

1. 8 

2.4 

1.8 

Average 

1. 1 

1.7 

2.2 



1. 6 

Table E-2 Variation of x^^ (x 10 

n 

° cm) with T and c 

' O 0 


970 

980 

990° C 

Average 

0. 01 

0. 5 

0.8 

1.3 

0. 9 

0. 02 

1. 2 

1. 9 

2. 6 

1. 9 

0.03 

1. 9 

3. 0 

4. 2 

3. 0 

0. 04 

2. 8 

4.3 

5. 8 

4.3 

Average 

1. 6 

2. 5 

3. 5 

2. 5 
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Multiple regression analyses (Refs. E4 and E5) of these data gave the following 
prediction equations; 

1) for the time in ms to indicate solidification, t^^^ in ms 




1.030 + 21. 4C ±0.60 
o 


= -59. 81 + 0. 0626 T ±0.34 

o 


t = -60. 34 + 21. 4 C + 0. 0626 T ± 0. 23 
bs o o 


0.692 + 55. 2 C - 675C^±0. 66 
o o 


^s 

= -1277. 18 + 2. 55 T - 0. 00127 
bs o o 


0.37 


6 2 

- 60. 68 + 55. 2 C. - 675 C_ + 0. 0626 ± 0. 25 


*bs 

‘bs 


o 


= -1278 + 21. 4 C +2.55T - 0. 00127 T ±0.25 
0 0 o 

= -1277. 06 + 55. 2 C - 676 C ^ + 2. 55 -0. 00127 T 

O O 0 c 

± 0.26 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

(7) 


x} = -0. 333 + 114. 3 C ±0. 92 
DS o 


XT = -89. 35 to 0. 0938 T ±1.45 

DS O 

3 

= 92. 21 + 114. 3 C t 0. 0938 T 

o o 

= 0. 04174 + 85. 17 C + 583. 2 C ^ ± 1. 02 


o 


o 


o 


X. ^ = 278. 12 -0. 656 T + 0. OOO: 
bs o 

X 6 

= -91. 92 + 85.17 C + 583. 2 
bs o 

xj = 273. 90 - 0. 6535 T + 0. 00 

DS O 

= 277.80 + 85.16 C + 583. 37 C ^ - 0. 6608 T + 0. 00038 
l 3 S O 0 o o 

± 0. 53 (16) 



(8) 

efore solidification, x, 

’ bs 

in 10“* 


(9) 


(10) 

± 0.42 

(11) 

± 1. 02 

(12) 

^ ± 1.60 

(13) 

+ 0. 0938 T ±0.47 
0 

(14) 

' ^ + 114.3 C ± 0.47 
0 o 

(15) 
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The first two prediction equations in each case relates t^g or x^s to one indepen- 
dent variable, Cq or Tq. The third equation is a linear prediction equation involving 
both Cq and Tq. The last five equations involve the squared Co terms, squared Tq 
terms, or both. The numbers after the "+"sign in each equation is the "adjusted 
standard error" (Refs. E4 and E5). 

Thus, tt,s can be computed from known values of Cq and Tq according to Eq. 3 
with an error of 0. 23. 90% of the time, the computed tj^g value of (Y estimated) is off 
from the true value of (Y observed) by 2 x 0.23 = 0.46 ms. as also can be seen from 
Table E-3. 

Similarly x^jg can be predicted from known values of and T , according to 
Eq. 11, with an error of 0. 42. 90% of the time, the computed x^^ value (Y estimated) 
is off from the true value (Y observed) by 2 x 0. 42 = 0 . 84 sec. , as can also be seen 
from Table E-3. 


Table E-3 Summary of Statistical Significances 


For t, 
bs 


Eq. 

Independent 

Variable 

Adjusted 

Standard 

Error 

Adjusted 

Regression 

Coefficient 

Variance 
Ratio, F 

Significance 

Level 

1 

C 

o 

0. 60 ms 

0.4014 

1. 921 

<90% 

2 

T 

o 

0.34 

0.8566 

27. 56 

99. 9% 

3 

C ,T 
o o 

0.23 

0. 9404 

38.31 

99. 9% 

4 

C ,C 2 
0 0 

0.66 

0.3022 

0.948 

<90% 

5 

T , T 2 
O 0 

0.37 

0. 8474 

13.06 

99% 

6 

C ,C 2,t 
o o o 

0. 25 

0. 9419 

26.21 

99. 9% 

7 

C T ,T 2 
0 o’ 0 

0.25 

0. 9401 

25.37 

99.9% 

8 

C ,C ^,T ,T 2 
0 0 0 o 

0.26 

0. 9417 

19.51 

99. 9% 









For 


Eq. 

Independent 

Variable 

Adjusted 

Standard 

Error 

Adjusted 

Regression 

Coefficient 

Variance 
Ratio, F 

Significance 

Level 

9 

C 

o 

0.92 xlO“® 

0.8353 

23.08 

99.9% 

10 

T 

o 

1.45 

0.5002 

3.34 

90% 

11 

C ,T 
o o 

0.42 

0. 9709 

81.94 

99.9% 

12 

C ,C 2 
o o 

1.02 

0.8180 

10.46 

95% 

13 

T ,T 2 

o o 

1.60 

0.4187 

1. 503 

<90% 

14 

C , C 2, T 
o o 0 

0.47 

0.9686 

50. 02 

99.9% 

15 


0.47 

0. 9678 

48. 69 

99. 9% 

16 

C ,C T ,T^ 

o’ 0 ’ 00 

0.53 

0. 9647 

32.92 

99. 9% 


Table E-3 gives the different sets of independent variables used for the above 
16 prediction equations for tj^^ and x, Also tabulated are the adjusted standard er- 
rors, adjusted regression coefficients, the values of F ratio, and the significance 
levels (Refs. E3 and E4). 

One sees again frona this table that the time to solidification t, depends not 
much on C , but strongly on T . The best prediction equation is a lifiear equation in 
C and T ,°i. e. , Eq. 3 (error^. 23 ms). On the other hand, the loss of material due 
to^ evaporation before solidification depends strongly on C , but not much on T , al- 
though the best prediction equation for x. is also a lineaf* equation in C and ® 

T^ i. e. , Eq. 11 (error 42 A). ® ° 

In designing or conducting space processing experiments, it is important to know 
the time of evaporation before solidification t, . In many cases, this time ma^ be too 
long to fit into, e. g. , a rocket program, or tne surface composition may be too much 
changed from the intended initial composition to be useful, hi other cases, this time 
t, may be too short for proper instrumentation and recording, e. g. , the beginning of 
sofidification or nucleation. Eq. 1-8 given above may, therefore, be used for Ge-Sb 
alloys within the initial concentration & temperature ranges studied. Of course, other 
ranges of the same or different alloys can also be similarly studied using the same 
models and techniques described above. 

The material evaporated away before solidification x^^ is also important in certain 
critical space processing operations such as the manufacture of precision parts includ- 
ing thin films. 
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